%\documentclass{Thesis}
\documentclass[a4paper, 12pt, oneside]{report}
%\documentclass{memoir}
\usepackage[utf8]{inputenc}
\usepackage{amsmath,amssymb,fancyhdr,url}
%\usepackage{amsfonts,amssymb,amsmath,graphicx,hyperref}
\usepackage[pdftex]{graphicx}
\usepackage{tikz}
\usetikzlibrary{shapes.geometric}
\usepackage{mathtools}
\usepackage{cite}
\bibliographystyle{plain}
\usepackage{lipsum}
\usepackage{eqnarray}
\usepackage[nottoc]{tocbibind}
\usepackage{setspace}
\usepackage{amsthm}
\usepackage[english]{babel}
\usepackage{rotating}
%\newtheorem{theorem}{Theorem}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}




\begin{document}


\begin{titlepage}
	\begin{center}
		\vspace*{1cm}
		
		\textbf{EVEN-ODD CYCLED HIGH-ORDER S-FDTD METHOD FOR MAXWELL'S EQUATIONS AND APPLICATION TO COPLANAR WAVEGUIDES}
		
%		\vspace{0.5cm}

		
		\vspace{2cm}
		
		\textbf{MANINDER SARAI}
		

		\vspace{3cm}
		
		A THESIS SUBMITTED TO\\
		THE FACULTY OF GRADUATE STUDIES\\
	 IN PARTIAL FULFILLMENT OF THE REQUIREMENTS\\
		FOR THE DEGREE OF \\
      MASTER OF SCIENCE
		
		\vspace{2cm}
		
		
		GRADUATE PROGRAM IN	Applied \& Industrial Mathematics\\
		YORK UNIVERSITY\\
		TORONTO, ONTARIO\\
		
		\vspace{1.5cm}
		
		April 2017\\
		
		\vspace{0.8cm}
		
		\textcopyright  Maninder Sarai, 2017
		
	\end{center}
\end{titlepage}
\pagenumbering{roman}
\newpage


\doublespacing

\begin{abstract}
	\addcontentsline{toc}{chapter}{Abstract}
	\thispagestyle{plain}
	\pagenumbering{roman}
	\setcounter{page}{2}
	In this thesis, a new even-odd cycled high-order splitting finite difference time domain scheme for Maxwell's equations in two dimensions
	is developed. The scheme uses fourth order spatial difference operators
	and even-odd time step technique to make it more accurate in both space and time. The scheme is energy-conserved, unconditionally stable and very efficient in computation. We analyze in detail the stability, dispersion and phase error for the scheme. We also prove its energy conservation, convergence of high order accuracy and convergence
	of divergence free approximation. Numerical experiments confirm the theoretical analysis results. Further, the developed scheme is applied
	to computations of the grounded coplanar waveguides, the elevated CPW
	and the complex transitions between CPW and rectangular waveguides.
\end{abstract}

\newpage

\begin{center}
	{\large \textbf{Acknowledgments}}
	\addcontentsline{toc}{chapter}{Acknowledgments}
	\thispagestyle{plain}
	\pagenumbering{roman}
	\setcounter{page}{3}
\end{center}
\normalsize	
\vspace{2cm}
	I acknowledge, with gratitude, Professor Dong Liang for guidance, supervision and encouragement throughout my years working on this thesis. Without him, this work would not have been completed. \\$ $\\
	
	I would like to thank my husband, parents for their continuous moral support and encouragement throughout my years of study and through the process of researching and writing this thesis.
	
	

\newpage

\tableofcontents
 \newpage
\listoftables
\newpage
 \listoffigures 
\newpage
\newpage


\newpage
\newpage
\pagenumbering{arabic}
\chapter{Introduction}

Maxwell's equations are widely used in computational electromagnetic applications such as, radio frequency, microwave, antennas, and air-craft radars and so on. Several ADI and splitting finite difference time domain methods have been developed to compute the solutions of Maxwell's equations. When computing modern problems of long distance wave propagations and moderately high frequency propagations, there are great interests to develop time and spatial high-order schemes.\\

Finite difference time domain method (FDTD) for Maxwell's equations was first introduced by Yee \cite{Yee1966} in 1966 which was later developed by other researchers \cite{Monk1994} \cite{Taflove2000} to a very efficient numerical algorithm in computational electromagnetics. However, the FDTD method is only conditionally stable and has large computational cost. Papers \cite{Namiki1999},  \cite{Zheng2000} proposed ADI-FDTD schemes for Maxwell's equations which are unconditionally stable and of second order accuracy. Further study was done in paper \cite{welfer2007}. Paper \cite{Chen2008} proposed EC-S-FDTDI and EC-S-FDTDII which are efficient and easy to compute. The schemes are energy conserved, unconditionally stable and non-dissipative. EC-S-FDTDI is first order in time and second order in space while EC-S-FDTDII is second order in both space and time. More recently, paper \cite{Liang2013} developed spatial fourth order energy conserved S-FDTD schemes, EC-S-FDTD(1,4) and EC-S-FDTD(2,4), for solving Maxwell's equations, which are fourth order accurate in space. In this thesis, we will develop a time second-order scheme for solving Maxwell's equations in two dimensions where the spatial accuracy is of fourth order. \\

In Chapter 2, we propose a new even-odd cycled splitting finite difference time domain scheme, Even-Odd cycled 4th order EC-S-FDTD, shorten to the EO-4th-EC scheme, with fourth order accuracy in space and second order accuracy in time. For EO-4th-EC scheme, we apply the spatial fourth order difference operators to a two stage splitting technique for each time step. The scheme consists of odd and even time steps where for the odd time step, electric field in y direction $[E_y]$ and the intermediate value of magnetic field $[H_z]$ are computed in Stage 1, following that electric field in x direction $[E_x]$ and $[H_z]$ are solved in Stage 2. For the even time step, $E_x$ and intermediate value $H_z$ are computed in Stage 1 and $E_y$ and $H_z$ are computed in Stage 2. The spatial fourth order difference operators are obtained by a linear combination of two central differences one with a spatial step and one with three spatial steps. The boundary node difference operators are carefully defined keeping in mind the energy conservation and fourth order accuracy in space. Another important feature is that the use of even-odd two cycles achieves high-order accuracy in time while only using two stages in EO-4th-EC. We analyze in detail the stability, dispersion and phase error for the scheme. We will prove the energy conservation, high order accuracy and divergence free approximation of the scheme. \\

Sections 2.1 and 2.2 lay out specifics of the model, proposed scheme, the initial and the boundary conditions. In Sections 2.3 and 2.4, the stability and dispersion analysis are analyzed in details. To find the stability of the scheme, the equivalent expressions for the even and odd time steps are computed by eliminating the intermediate terms. The expressions further allow us to compute the growth factor for each time step and the scheme overall.  The growth factors help to determine the dispersion relationships of our scheme at each time step. Sections 2.5 and 2.6 comprise of theoretical analysis of energy conservation and convergence. Section 2.7 focuses on numerical dispersion results and numerical experiments. The growth factors are used to compare the phase velocity error of the proposed scheme to other schemes such as ADI-FDTD, CN, EC-S-FDTDI, EC-S-FDTDII, EC-S-FDTD (2,4). The energy conservation, accuracy errors and the divergence free approximations are numerically computed and compared.\\

Chapter 3 combines experimental results of different domains along with various types of sources for the electromagnetic wave propagations. Section 3.1 lays out the model and scheme inclusive of a source term. Section 3.2 looks at the behaviour of initial magnetic/electric field vs a magnetic/electric current source. Different domains are evaluated. Vacuum domain is compared with Perfectly electric conducting (PEC) boundary and Perfectly Matched layer(PML) and with the presence of a source/initial current. PEC is an idealized material exhibiting infinite electric conductivity and zero electric resistance. PML is made of lossy material which is used to absorb the fields traveling away from the interior of the grid. PML was first introduced by J.P. Berenger in 1994 \cite{Berenger1994}. In each domain, metallic layers are added to see the behaviour of the moving current in different layers. The source of electric/magnetic current is defined in a shape of a thin conductor.\\

Chapter 4 looks at the applications of the EO-4th-EC scheme to the grounded coplanar waveguides in complex structures of printed circuit boards. Coplanar waveguide is a type of transmission line which can be fabricated using printed circuit board type of technology. Coplanar waveguides are made of two parallel plates made of conducting material that run for a finite difference with some dielectric medium in the center. Coplanar waveguides are often compared to microstrip, rectangular waveguides, and striplines. However, the geometry of CPW results in minimum dielectric loss which is why it is considered to be the best transmission line. Although, in the past, some believed that CPW has a higher conductive loss. With more research, CPW resulted in the same conductor loss as its comparable microstrip and much lower dielectric loss. Also, it was found that conductor loss depends on the roughness of the surface rather than the geometry. Coplanar Waveguides are cheaper and less bulkier to make to use for printed circuit boards. The fabrication of CPW is often tailored to get the optimal result. The waves inside Coplanar waveguides normally travel at a very high frequency. Electromagnetic waves travel very efficiently due to CPW's geometry ending up in low losses such as dielectric, conductor and radiation loss.  Experimental results are run in MATLAB to confirm the behaviour of waves in different sets of domain. \\

Section 4.1 lists the details of many applications in regards to a coplanar waveguide. Section 4.2 looks at one of many ways to excite a coplanar waveguide. In the section, we look at the behaviour of EM waves when there are two lumped ports carrying current towards or away from the center conductor. Section 4.3 focuses more on a grounded coplanar waveguide when there are two sources inside in a shape of a cylindrical wire carrying current. Current moving in the same direction and the opposite direction are compared. In Section 4.4, an elevated coplanar waveguide is simulated as it is known to have a wider bandwidth for the EM waves to move as seen in paper \cite{Hofschen1996}. Section 4.4 focuses on a transition between a coplanar waveguide and a rectangular waveguide from \cite{Deslandes2001}. This chapter also compares the  dielectric and conductor loss inside a coplanar waveguide and a microstrip. It also compares the dielectric and conductor loss inside a CPW and an elevated CPW. \\

 Overall, the proposed scheme is found to be unconditionally stable and non-dissipative. The scheme also conserves energy, has higher accuracy than previous schemes. This high order scheme can be used to solve applications of MMIC such as coplanar, rectangular waveguides and dipole antennas. Coplanar and rectangular waveguides are commonly used in electrical printed circuit boards. We look deeply into the systems of coplanar waveguides and how the waves travel in different media. Grounded coplanar waveguide, elevated coplanar waveguide and transition between different waveguides are further analyzed.
 

\newpage
\singlespacing
\chapter{The EO-EC-S-FDTD Scheme and Theoretical Analysis}
\doublespacing
In this chapter, we propose the even-odd 4th order energy conserved splitting finite difference time domain scheme for Maxwell's equations. The chapter first provides details on stability and dispersion analysis of the scheme. Theoretical analysis of energy conservation and convergence is then considered in this chapter. Numerical results are presented in the last section and they show that the scheme is of second order in time and fourth order in space and is energy conserved.
\section{Maxwell's Equations in 2D}
Consider the 2D transverse electric polarization case in a lossless medium and there is no source. We have\\
\begin{equation}
 \frac{\partial{E_x}}{\partial{t}}= \frac{1}{\epsilon} \frac{\partial{H_z}}{\partial{y}}
\end{equation}
\begin{equation}
 \frac{\partial{E_y}}{\partial{t}}= -\frac{1}{\epsilon} \frac{\partial{H_z}}{\partial{x}} \\
\end{equation}
\begin{equation}
\frac{\partial{H_z}}{\partial{t}}= \frac{1}{\mu}(\frac{\partial{E_x}}{\partial{y}}-\frac{\partial{E_y}}{\partial{x}} ),
\end{equation}
where, $\epsilon$ and $\mu$ are the permittivity and permeability, $E=(E_x(x,y,t)), E_y(x,y,t))$ is the electric field and $H_z=H_z(x,y,t)$ is the magnetic field for $(x,y) \in \Omega$ $ =(0,a)\times (0,b)$ and $t \in (0,T]$.
Perfectly conducting boundary condition is satisfied on the boundary:
$(E,0) \times (\vec{n},0)$ = 0 on (0,T] $\times \partial{\Omega},$
 where $\partial{\Omega}$ is the boundary of $\Omega$ and $\vec{n}$ is the outward normal vector on $\partial{\Omega}$. The initial conditions are:
$\mathbf{E}(x,y,0)= \mathbf{E_0}(x,y)$=$(E_{x_0} (x,y),E_{y_0} (x,y))$, $H_z (x,y,0)=H_{z_0}(x,y).$\\

Let the partition of space domain $\Omega$ and the time interval [0,T] be uniformly staggered grid:
\begin{multline}
x_i=i\Delta x, x_{i+\frac{1}{2}}=x_i +\frac{1}{2}\Delta x, i=0,1,..,I-1, x_I=I\Delta x, \nonumber \\
y_j=j\Delta y, y_{j+\frac{1}{2}}=y_j +\frac{1}{2}\Delta y, j=0,1,..,J-1, y_J=J\Delta y, \nonumber \\
t^n = n\Delta t, t^{n+\frac{1}{2}}=t^n+\frac{1}{2}\Delta t, n=0,1,..,N-1, t^N=N \Delta t \nonumber
\end{multline}
where $\Delta x$ and $\Delta y$ are the mesh sizes along the x and y directions, respectively, and $\Delta t$ is the time step size, I, J and N are integers. For a function F(t,x,y), we let $F^m_{\alpha,\beta}=F(m\Delta t, \alpha \Delta x, \beta\Delta y),$ and define\\
\begin{center}
\begin{eqnarray}
&& \delta_t F^m_{\alpha,\beta}=\frac{F^{m+\frac{1}{2}}_{\alpha,\beta}-F^{m-\frac{1}{2}}_{\alpha,\beta}}{\Delta t}, \qquad \delta_x F^m_{\alpha,\beta}= \frac{F^{m}_{\alpha+\frac{1}{2},\beta}-F^{m}_{\alpha-\frac{1}{2},\beta}}{\Delta x}, \nonumber \\
&& \delta_y F^m_{\alpha,\beta}=\frac{F^{m}_{\alpha,\beta+\frac{1}{2}}-F^{m}_{\alpha,\beta-\frac{1}{2}}}{\Delta y} \qquad \delta_u \delta_v F^m_{\alpha,\beta}=\delta_u(\delta_v F^m_{\alpha,\beta}),\nonumber \\ && u,v=x,y. \\ \nonumber
\end{eqnarray} 
\end{center}
For the strict interior nodes,we further define the fourth order difference operator to $\frac{\partial{E_y}}{\partial{x}}(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}})$ by
\begin{center}
$\Lambda_x E^n_{y_{i+\frac{1}{2},j+\frac{1}{2}}} = \frac{1}{8}(9\delta_x-\delta_{2,x})E^n_{y_{i+\frac{1}{2},j+\frac{1}{2}}},$\\
\end{center}
where,
\begin{center}
$\delta_{2,x} E^n_{y_{\alpha,\beta}}= \frac{E^{n}_{y_{\alpha+\frac{3}{2},\beta}}-E^{n}_{y_{\alpha-\frac{3}{2},\beta}}}{3\Delta x}$ \\
\end{center}

However, $\delta_{2,x}E^{n}_{y_{\frac{1}{2}, j+\frac{1}{2}}}$ and $\delta_{2,x}E^{n}_{y_{I-\frac{1}{2}, j+\frac{1}{2}}}$ go out of domain. We can define


$$\tilde{\delta}_{2,x}E^n_{y_{\frac{1}{2},j+\frac{1}{2}}} = \frac{E^n_{y_{1,j+\frac{1}{2}}}+E^n_{y_{2,j+\frac{1}{2}}}-2E^n_{y_{0,j+\frac{1}{2}}}}{3\Delta x}$$ 
$$\tilde{\delta}_{2,x}E^n_{y_{I-\frac{1}{2},j+\frac{1}{2}}} = \frac{2E^n_{y_{I,j+\frac{1}{2}}}-E^n_{y_{I-1,j+\frac{1}{2}}}-E^n_{y_{I-2,j+\frac{1}{2}}}}{3\Delta x}$$


where the following relations in [3] are used on the boundary nodes:
$$E_y (x_{-1},y_{j+\frac{1}{2}},t)=2 E_y (x_{0},y_{j+\frac{1}{2}},t)-E_y (x_{1},y_{j+\frac{1}{2}},t)+O(\Delta x^5)$$
$$E_y (x_{I+1},y_{j+\frac{1}{2}},t)=2 E_y (x_{I},y_{j+\frac{1}{2}},t)-E_y (x_{I-1},y_{j+\frac{1}{2}},t)+O(\Delta x^5)$$

For the near boundary points with i=1 and i=I-1,
$$\tilde{\Lambda}_x E^n_{y_{\frac{1}{2},j+\frac{1}{2}}} = \frac{1}{8}(9\delta_x-\tilde{\delta}_{2,x})E^n_{y_{\frac{1}{2},j+\frac{1}{2}}},$$
$$\tilde{\Lambda}_x E^n_{y_{I-\frac{1}{2},j+\frac{1}{2}}} = \frac{1}{8}(9\delta_x-\tilde{\delta}_{2,x})E^n_{y_{I-\frac{1}{2},j+\frac{1}{2}}}$$

for j= 0,1,...,J-1.


Similarly, we can define $\Lambda_y F^m_{\alpha,\beta}$ in the y-direction.
\\

We can split the Maxwell's equations in each time interval $[t^n, t^{n+1}]$ as following:
\begin{equation}
\frac{\partial{E_x}}{\partial{t}}=\frac{1}{\epsilon}\frac{\partial{H_z}}{\partial{y}},\\
\qquad    \frac{1}{2}\frac{\partial{H_z}}{\partial{t}}=\frac{1}{\mu}\frac{\partial{E_x}}{\partial{y}} \end{equation}
and 
\begin{equation}
\frac{\partial{E_y}}{\partial{t}}=-\frac{1}{\epsilon}\frac{\partial{H_z}}{\partial{x}},\\
\qquad \frac{1}{2}\frac{\partial{H_z}}{\partial{t}}=-\frac{1}{\mu}\frac{\partial{E_y}}{\partial{x}}
 \end{equation}


\section{The EO-4th-EC Scheme}
We propose the Even-Odd cycled 4th Order EC-S-FDTD scheme, shorten to EO-4th-EC. The proposed scheme is a combination of EC-S-FDTD(1,4) and the Even-Odd cycle technique. \\

\textbf{Stage 1: The Odd time step.}\\
\textbf{At every odd time step, from $t^{(2k)}$ to $ t^{(2k+1)}$, use $E^{(2k)}_x$, $E^{(2k)}_y$ and $H^{(2k)}_z$ to compute $E^{(2k+1)}_x$, $E^{(2k+1)}_y$ and $H^{(2k+1)}_z$ by the fourth-order spatial scheme.}\\
\textit{Stage 1.1: First, compute $E^{(2k+1)}_y$ and intermediate variable $H^*_z$ using $H^{(2k)}_z$ and $E^{(2k)}_y$; for j=0,1,..,J-1}\\ 
For the interior points, i=2,3,..,I-2
\begin{multline}
\frac{E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}-E^{(2k)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t}
= - \frac{1}{2\epsilon}\varLambda_x\{H^*_{z_{i,j+\frac{1}{2}}} + H^{(2k)}_{z_{i,j+\frac{1}{2}}}\}, \\
 \frac{H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t}
= - \frac{1}{2\mu}\varLambda_x\{E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\},
\end{multline}
For boundary nodes, i=1 and i=I-1,
\begin{multline}
\frac{E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}-E^{(2k)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t}
= - \frac{1}{2\epsilon}\widetilde{\varLambda_x}\{H^*_{z_{i,j+\frac{1}{2}}} + H^{(2k)}_{z_{i,j+\frac{1}{2}}}\}, \\
\frac{H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t}
 = - \frac{1}{2\mu}\widetilde{\varLambda_x}
\{E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\},
\end{multline}
\textit{Stage 1.2: Then compute $E^{(2k+1)}_x$ and $H^{(2k+1)}_z$ using $H^{(*)}_z$ and $E^{(2k)}_x$, for i = 0,1,..,I-1}\\ 
For interior nodes, j=2,3,..,J-2,
\begin{multline}
\frac{E^{(2k+1)}_{x_{i+\frac{1}{2},j}}-E^{(2k)}_{x_{i+\frac{1}{2},j}}}{\varDelta t}
 =  \frac{1}{2\epsilon}\varLambda_y\{H^{(2k+1)}_{z_{i+\frac{1}{2},j}} + H^{*}_{z_{i+\frac{1}{2},j}}\}, \\
\frac{H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t}  =\frac{1}{2\mu}\varLambda_y\{E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\},
\end{multline}
For boundary nodes,j=1 and j=J-1,
\begin{multline}
\frac{E^{(2k+1)}_{x_{i+\frac{1}{2},j}}-E^{(2k)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} = \frac{1}{2\epsilon}\widetilde{\varLambda_y}\{H^{(2k+1)}_{z_{i+\frac{1}{2},j}} + H^{*}_{z_{i+\frac{1}{2},j}}\},\\
\frac{H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu}\widetilde{\varLambda_y}\{E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\},
\end{multline}

\textbf{Stage 2. The Even time step.}\\

\textbf{At every even time step, from $t^{(2k+ 1)}$ to $t^{(2k+2)}$, use $E^{(2k+1)}_x$, $E^{(2k+1)}_y$ and $H^{(2k+1)}_z$ to compute $E^{(2k+2)}_x$, $E^{(2k+2)}_y$ and $H^{(2k+2)}_z$ by firstly fourth order difference y-direction scheme and secondly fourth order difference x-direction scheme .} \\

\textit{Stage 2.1: First compute $E^{(2k+2)}_x$ and intermediate variable $H^{**}_z$ using $H^{(2k+1)}_z$ and $E^{(2k+1)}_x$, i=0,1,..,I-1.} \\
For interior nodes, j = 2,3,..,J-2
\begin{multline}
\frac{E^{(2k+2)}_{x_{i+\frac{1}{2},j}}-E^{(2k+1)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} =  \frac{1}{2\epsilon}\varLambda_y\{H^{**}_{z_{i+\frac{1}{2},j}} + H^{(2k+1)}_{z_{i+\frac{1}{2},j}}\}, \\
\frac{H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu}\varLambda_y\{E^{(2k+2)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\},
\end{multline}
For boundary nodes, j=1 and J-1.
\begin{multline}
\frac{E^{(2k+2)}_{x_{i+\frac{1}{2},j}}-E^{(2k+1)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} =  \frac{1}{2\epsilon}\widetilde{\varLambda_y}\{H^{**}_{z_{i+\frac{1}{2},j}} + H^{(2k+1)}_{z_{i+\frac{1}{2},j}}\},\\
\frac{H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu}\widetilde{\varLambda_y}\{E^{(2k+2)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\},
\end{multline} \\$ $\\

\textit{Stage 2.2: Then compute $E^{(2k+2)}_y$ and $H^{(2k+2)}_z$ using $H^{**}_z$ and $E^{(2k+1)}_y$, j=0,1,..,J-1}\\
For interior nodes, i=2,3,..,I-2.
\begin{multline}
\frac{E^{(2k+2)}_{y_{i,j+\frac{1}{2}}}-E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t} =  -\frac{1}{2\epsilon}\varLambda_x\{H^{(2k+2)}_{z_{i,j+\frac{1}{2}}} + H^{**}_{z_{i,j+\frac{1}{2}}}\}, \\
\frac{H^{(2k+2)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} =-\frac{1}{2\mu}\varLambda_x\{E^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\},
\end{multline}
For boundary nodes, i=1, and i = I-1:
\begin{multline}
\frac{E^{(2k+2)}_{y_{i,j+\frac{1}{2}}}-E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t} =  -\frac{1}{2\epsilon}\widetilde{\varLambda_x}\{H^{(2k+2)}_{z_{i,j+\frac{1}{2}}} + H^{**}_{z_{i,j+\frac{1}{2}}}\},\\
\frac{H^{(2k+2)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} =-\frac{1}{2\mu}\widetilde{\varLambda_x}\{E^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\},
\end{multline}

The boundary conditions are given by
$$E^*_{x_{i+\frac{1}{2},0}}=E^*_{x_{i+\frac{1}{2},J}}=E^{(2k+2)}_{x_{i+\frac{1}{2},0}}=E^{(2k+2)}_{x_{i+\frac{1}{2},J}}=E^{(2k+2)}_{y_{0,j+\frac{1}{2}}}=E^{(n+1)}_{y_{I,j+\frac{1}{2}}} = 0,$$
and the initial conditions are given by
\begin{center}
$E^0_{x_{\alpha,\beta}}= E_{x0}(\alpha \Delta x, \beta \Delta y); E^0_{y_{\alpha,\beta}} = E_{y0}(\alpha \Delta x, \beta \Delta y);$ \\$ H^0_{z_{\alpha,\beta}}= H_{z0}(\alpha \Delta x, \beta \Delta y).$
\end{center}

\newpage
\section{Dispersion and Stability Analysis}

We study the stability analysis and the dispersion analysis for our EO-4th-EC scheme. According to the Fourier Analysis, we define that for the time step (2k) to (2k+1),
\begin{multline}
E^{(2k+1)}_{\alpha,\beta}=E^{(2k)} \xi^{(1)} e^{-i(k_x\alpha \Delta x + k_y\beta\Delta y)},\\
H^{(2k+1)}_{\alpha,\beta}=H^{(2k)}_z \xi^{(1)} e^{-i(k_x\alpha \Delta x + k_y\beta\Delta y)}
\end{multline}
And for time step (2k+1) to (2k+2), we let:
\begin{multline}
E^{(2k+2)}_{\alpha,\beta}=E^{(2k+1)}\xi^{(2)} e^{-i(k_x\alpha \Delta x + k_y\beta\Delta y)},\\
H^{(2k+2)}_{\alpha,\beta}=H^{(2k+1)}_z \xi^{(2)} e^{-i(k_x\alpha \Delta x + k_y\beta\Delta y)}.
\end{multline}
Here, $\xi^{(1)}$ and $\xi^{(2)}$ are the growth factors for each time step. Growth factors are complex time eigenvalues whose magnitude will determine the stability and dissipation properties of the numerical scheme. \textbf{K}=$(k_x,k_y)$ is the real wave number of the arbitrary harmonic wave component.\\


Equivalent expressions are formed by eliminating the intermediate terms and have them in simplified form. They are further used to compute the numerics such as the stability, dispersion and phase error. We will determine the equivalent expression for each time separately.
\subsubsection{Equivalent Expression for the EO-4th-EC Scheme}\\
For the odd time step, it gets\\
\begin{equation} \frac{E^{(2k+1)}_{x_{i+\frac{1}{2},j}}-E^{(2k)}_{x_{i+\frac{1}{2},j}}}{\varDelta t}  =  \frac{1}{2\epsilon}\varLambda_y\left\{H^{(2k+1)}_{z_{i+\frac{1}{2},j}} + H^{(2k)}_{z_{i+\frac{1}{2},j}}\right\} -\frac{\Delta t}{4\mu\epsilon}{\varLambda_x}{\varLambda_y}\left\{E^{(2k+1)}_{y_{i+\frac{1}{2},j}}+E^{(2k)}_{y_{i+\frac{1}{2},j}}\right\},
\end{equation}
\begin{equation}
 \frac{E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}-E^{(2k)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t}  
 = -\frac{1}{2\epsilon}{\varLambda_x}\left\{H^{(2k+1)}_{z_{i,j+\frac{1}{2}}} + H^{(2k)}_{z_{i,j+\frac{1}{2}}}\right\} 
 +\frac{\Delta t}{4\mu\epsilon}\varLambda_x\varLambda_y \left\{E^{(2k+1)}_{x_{i,j+\frac{1}{2}}}+E^{(2k)}_{x_{i,j+\frac{1}{2}}}\right\},
 \end{equation}
 \begin{equation}
 \frac{H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t}  =   \frac{1}{2\mu} \left\{\varLambda_y \left\{E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\right\} \right. -\varLambda_x \left. \left\{ E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} \right\}\right\} .
\end{equation}

For the even time step, it gets
\begin{equation}
 \frac{E^{(2k+2)}_{x_{i+\frac{1}{2},j}}-E^{(2k+1)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} =  \frac{1}{2\epsilon}{\varLambda_y}\left\{H^{(2k+2)}_{z_{i+\frac{1}{2},j}} + H^{(2k+1)}_{z_{i+\frac{1}{2},j}}\right\}  +\frac{\Delta t}{4\mu\epsilon}{\varLambda_x}{\varLambda_y}\left\{E^{(2k+2)}_{y_{i+\frac{1}{2},j}}+E^{(2k+1)}_{y_{i+\frac{1}{2},j}}\right\}
 \end{equation}
 \begin{equation}
  \frac{E^{(2k+2)}_{y_{i,j+\frac{1}{2}}}-E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t} =  -\frac{1}{2\epsilon}{\varLambda_x}\left\{H^{(2k+2)}_{z_{i,j+\frac{1}{2}}} + H^{(2k+1)}_{z_{i,j+\frac{1}{2}}}\right\} -\frac{\Delta t}{4\mu\epsilon}\varLambda_x\varLambda_y\left\{E^{(2k+2)}_{x_{i,j+\frac{1}{2}}}+E^{(2k+1)}_{x_{i,j+\frac{1}{2}}}\right\}
 \end{equation}
 \begin{equation}
 \frac{H^{(2k+2)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu}\{\varLambda_y\left\{E^{(2k+2)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\right\}-\varLambda_x\left\{ E^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\}\right\}.
 \end{equation}

For the odd time step, we substitute equation (2.15) into equivalent expressions (2.17)-(2.19) to analyze the stability for which we get the following system:
 $$AX^{(2k)}=0,$$ Where, $X^{(2k)} = [E^{(2k)}_x, E^{(2k)}_y, H^{(2k)}_z]^T$, given in (2.15), and A is the coefficient matrix (2.23) for the odd time step with\\
 $$ u_x=\frac{\sin{(\frac{3k_x \Delta x}{2})}-27\sin{(\frac{k_x \Delta x}{2})} }{24\Delta x},$$
  $$v_y=\frac{\sin{(\frac{3k_y \Delta y}{2})}-27\sin{(\frac{k_y \Delta y}{2})} }{24\Delta y}.$$

\begin{table*}
	\begin{equation}
	A= \begin{bmatrix}
	\xi^{(1)} - 1 & -\frac{\Delta t^2}{\mu\epsilon}u_x v_y (\xi^{(1)} +1)& -i\frac{\Delta t}{\epsilon}v_y(\xi^{(1)} +1)\\
	\frac{\Delta t^2}{\mu\epsilon}u_xv_y(\xi^{(1)} +1) & (\xi^{(1)} -1) & i \frac{\Delta t}{\epsilon}u_x(\xi^{(1)} +1)\\
	-i\frac{\Delta t}{\mu}v_y(\xi^{(1)} +1) & i \frac{\Delta t}{\mu}u_x(\xi^{(1)} +1) & (\xi^{(1)} -1)
	\end{bmatrix}
	\end{equation}
\end{table*}

 In order to find the growth matrix of the odd time step, we solve det(A)=0. The following equation is obtained:
 \begin{equation}
 (\xi^{(1)} -1)[{Q \xi^{(1)}}^{2} +2R\xi^{(1)}+Q]=0
 \end{equation}
 where,
 \begin{equation}
 Q=1+\frac{\Delta t^2}{\mu\epsilon}(u^2_x+v^2_y)+\frac{\Delta t^4}{\mu^2\epsilon^2}u^2_x v^2_y,
 \end{equation}
 \begin{equation}
 R=-1+\frac{\Delta t^2}{\mu\epsilon}(u^2_x+v^2_y)+\frac{\Delta t^4}{\mu^2\epsilon^2}u^2_x v^2_y.
 \end{equation}

The roots of (22) are $\xi^{(1)}_1 = 1$ and $\xi^{(1)}_{2,3}=-\frac{R}{Q}\pm i\frac{ \sqrt{Q^2 -R^2}}{Q} \rightarrow |\xi^{(1)}_{2,3}|=1$ \\
$~$
\\
For the even time step, similarly substituting equation (2.16) into equivalent expressions (2.20)-(2.22), we get
$$BX^{(2k+1)}=0$$
Where, $X^{(2k+1)} = [E^{(2k+1)}_x, E^{(2k+1)}_y, H^{(2k+1)}_z]^T$, given in (2.16), and B is the coefficient matrix (2.27) for the even time step.
\begin{table*}
	\begin{equation}
	B= \begin{bmatrix}
	\xi^{(2)} - 1 & \frac{\Delta t^2}{\mu\epsilon}u_x v_y (\xi^{(2)} +1)& -i\frac{\Delta t}{\epsilon}v_y(\xi^{(2)} +1)\\
	-\frac{\Delta t^2}{\mu\epsilon}u_x v_y(\xi^{(2)} +1) & (\xi^{(2)} -1) & i \frac{\Delta t}{\epsilon}u_x(\xi^{(2)} +1)\\
	-i\frac{\Delta t}{\mu}v_y(\xi^{(2)} +1) & i \frac{\Delta t}{\mu}u_x(\xi^{(2)} +1) & (\xi^{(2)} -1)
	\end{bmatrix}
	\end{equation}
\end{table*}
Solving for det(B)=0, we get the following equation:
\begin{equation}
(\xi^{(2)} -1)[{Q \xi^{(2)}}^{2} +2R\xi^{(2)}+Q]=0
\end{equation}
Here, the roots are $\xi^{(2)}_1=1, \xi^{(2)}_{2,3} = -\frac{R}{Q}\pm i\frac{ \sqrt{Q^2 -R^2}}{Q} \rightarrow |\xi^{(2)}_{2,3}|=1$ \\
$~$
\\
The growth factor for the scheme combining Stage 1 with Stage 2 is $\xi= \xi^{(1)} \xi^{(2)}$. By calculating, we have that\\
\begin{equation}
   \xi = \left\{
     \begin{array}{lr}
      \frac{2R^2-Q^2}{Q^2} + i \frac{2R \sqrt{Q^2-R^2}}{Q^2} \\
       1 \\
       \frac{2R^2-Q^2}{Q^2} - i \frac{2R \sqrt{Q^2-R^2}}{Q^2}.
     \end{array}
   \right.
\end{equation}

It is clear that $|\xi| =1$.\\
Thus we have the following theorem. 

\begin{theorem}[]
The EO-4th-EC Scheme is non-dissipative and unconditionally stable.
\end{theorem}


\section{Dispersion Relationship}
The solution $\xi$ of the characteristic polynomial equation (2.25) and (2.28) determines the dispersion properties of the scheme. $|\xi|$ determines the amplitude (or dissipative) error, for which arctan(Im$\xi$/Re$\xi$) determines phase (or dispersive error) which is found in Section 2.7 .

\subsubsection{Dispersion Relationship for the EO-4th-EC Scheme}
First, we can find the dispersion relationship for the EO-4th-EC scheme. We let $\gamma=j \omega$ and let
\begin{equation}
\xi^{(1)}= \xi^{(2)}= \xi_{l}= e^{\gamma\Delta tn}
\end{equation}
Since the growth factors for the odd and even time step are same, we get the same dispersion relationships for both even and odd time steps. To find the dispersion relationship for each time step, we consider the part with complex roots of the stability equations (2.25) and (2.28):
\begin{equation}
{Q \xi_{l}}^{2} +2R\xi_{l}+Q=0
\end{equation}
By simplifying and dividing every term by $\xi^{(l)}$ and we get:
\begin{equation}
Q( \xi_{l}+\xi_{l}^{-1}) +2R=0 %needs to be fixed
\end{equation}

Using the definition (2.30), we get:
\begin{equation}
\xi_{l}+\xi_{l}^{-1}=-4\sin^2(\frac{\omega\Delta t}{2})+2
\end{equation}

Substituting equation (2.33), (2.26) and (2.27) into equation (2.32), we get,
\begin{eqnarray}
&& [1+ \frac{\Delta t^2}{\mu\epsilon}(u^2_x + v^2_y)+\frac{\Delta t^4}{\mu ^2\epsilon ^2}u^2_x v^2_y](-4\sin^2(\frac{\omega\Delta t}{2})+2)\nonumber\\
&& \quad +2[-1+ \frac{\Delta t^2}{\mu\epsilon}(u^2_x + v^2_y)+\frac{\Delta t^4}{\mu ^2\epsilon ^2}u^2_x v^2_y]=0
\end{eqnarray}

Simplifying we get the dispersion relationship,
\begin{equation}
\left[\frac{\Delta t^2}{\mu\epsilon}(u^2_x + v^2_y)+\frac{\Delta t^4}{\mu ^2\epsilon ^2}u^2_x v^2_y\right]\cos^2(\frac{\omega\Delta t}{2})=\sin^2(\frac{\omega\Delta t}{2})
\end{equation}
Using the definition of $u_x$ and $v_y$,
\begin{eqnarray}
&& \left\{\frac{\Delta t^2}{\mu\epsilon}\left\{{\left(\frac{\sin{(\frac{3k_x \Delta x}{2})}-27\sin{(\frac{k_x \Delta x}{2})}}{24\Delta x}\right)}^2 \right. \right. \left. + {\left(\frac{\sin{(\frac{3k_y \Delta y}{2})}-27\sin{(\frac{k_y \Delta y}{2})}}{24\Delta y}\right)}^2\right\} \nonumber \\
&& \quad + \frac{\Delta t^4}{\mu ^2\epsilon ^2}{\left(\frac{\sin{(\frac{3k_x \Delta x}{2})}-27\sin{(\frac{k_x \Delta x}{2})} }{24\Delta x}\right)}^2 \left. {\left(\frac{\sin{(\frac{3k_y \Delta y}{2})}-27\sin{(\frac{k_y \Delta y}{2})} }{24\Delta y}\right)}^2\right\}\cos^2(\frac{\omega\Delta t}{2}) \nonumber \\
&& \qquad =\sin^2(\frac{\omega\Delta t}{2})
\end{eqnarray}

which can be simplified to the dispersion relationship for the EO-4th-EC Scheme as
\begin{eqnarray}
&& \left\{\Delta t^2 c^2\left\{{\left(\frac{\sin{(\frac{3k_x \Delta x}{2})}-27\sin{(\frac{k_x \Delta x}{2})}}{24\Delta x}\right)}^2 \right. \right. \left. + {\left(\frac{\sin{(\frac{3k_y \Delta y}{2})}-27\sin{(\frac{k_y \Delta y}{2})} }{24\Delta y}\right)}^2 \right\}+ \nonumber \\
&& \quad \Delta t^4 c^4  {\left(\frac{\sin{(\frac{3k_x \Delta x}{2})}-27\sin{(\frac{k_x \Delta x}{2})} }{24\Delta x}\right)}^2  \left. {\left(\frac{\sin{(\frac{3k_y \Delta y}{2})}-27\sin{(\frac{k_y \Delta y}{2})} }{24\Delta y}\right)}^2\right\} \cos^2(\frac{\omega\Delta t}{2}) \nonumber \\
&& \qquad =\sin^2(\frac{\omega\Delta t}{2})
\end{eqnarray}

where $c = \frac{1}{\sqrt{\mu\epsilon}}$.
Note, here $k_x = k\cos(\phi)$ and $k_y = k\sin(\phi)$. As $\Delta t, \Delta x, \Delta y \rightarrow 0$, the dispersion relationship converges to
$$k^2_x + k^2_y = \frac{\omega^2}{c^2}.$$



\section{Energy Conservation}

For grid functions defined on the staggered grid:\\ $ $ \\
$U:={{U_{i+\frac{1}{2},j}}}, V:={{V_{i,j+\frac{1}{2}}}}, W={{W_{i+\frac{1}{2},j+\frac{1}{2}}}}, 
\vec{F}:={({U_{i+\frac{1}{2},j}},V_{i,j+\frac{1}{2}})} $,\\ $ $\\
the discrete $L^2$ energy norms are used:\\ $ $\\
$\| U\|^2 _{E_x}= \sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} \mid U_{i+\frac{1}{2},j}\mid^2 \Delta x\Delta y,$ 
$\| V\|^2 _{E_y}= \sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} \mid V_{i+\frac{1}{2},j}\mid^2 \Delta x\Delta y$ \\ $ $\\
$\| W\|^2 _{H}= \sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} \mid W_{i+\frac{1}{2},j+\frac{1}{2}}\mid^2 \Delta x\Delta y$
$\| \vec{F}\|^2 _{E}= \| U\|^2 _{E_x}+ \| V\|^2 _{E_y} $\\

\begin{theorem}[Discrete Energy Conservations]

	\textit{For integer $n\ge0$, if $\vec{E^n}:={(E^n_{x_{i+\frac{1}{2},j}},E^n_{y_{i,j+\frac{1}{2}}})}$ and $H_z^n:={H^n_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}$ are the solutions of EO-4th-EC scheme, then the following discrete energies are constant.}
	\begin{equation}
	\|\epsilon^{\frac{1}{2}} \vec{\textbf{E}}^{n+1} \|^2_E +\|\mu^{\frac{1}{2}} {H_z}^{n+1} \|^2_H = \|\epsilon^{\frac{1}{2}} \vec{\textbf{E}}^{n} \|^2_E +\|\mu^{\frac{1}{2}} {H_z}^{n} \|^2_H
	\end{equation}
	where n=2k and 2k+1,
	and 
	\begin{multline}
	\|\epsilon^{\frac{1}{2}}  \frac{(\vec{\textbf{E}}^{n+2} - \vec{\textbf{E}}^{n})}{ (2\Delta t)} \|^2_E 
	+\|\mu^{\frac{1}{2}} \frac{({H_z}^{n+2} -{H_z}^{n})}{(2 \Delta t)} \|^2_H \\ =
	\|\epsilon^{\frac{1}{2}}  \frac{(\vec{\textbf{E}}^{n+1} - \vec{\textbf{E}}^{n-1})}{ (2\Delta t)} \|^2_E 
	+\|\mu^{\frac{1}{2}} \frac{({H_z}^{n+1} -{H_z}^{n-1})}{(2 \Delta t)} \|^2_H 
	\end{multline}
	where n=2k and 2k+1.
\end{theorem} 
\begin{proof}
Multiplying both sides of (2.7a) with $\epsilon \Delta t(E^{{2k+1}}_{i,j+\frac{1}{2}}+E^{{2k}}_{i,j+\frac{1}{2}})$ and Multiplying both sides of (2.7b) with $\mu \Delta t (H^{*}_{i+\frac{1}{2},j+\frac{1}{2}}+H^{2k}_{i+\frac{1}{2},j+\frac{1}{2}}) $, we can get:
\begin{equation}
\epsilon \left[ (E_y_{i,j+\frac{1}{2}}^{(2k+1)})^2 - (E_y_{i,j+\frac{1}{2}}^{(2k)})^2 \right] = -\frac{\Delta t}{2}\varLambda_x {(H^*_{z_{i+\frac{1}{2},j+\frac{1}{2}}} + H^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}})} (E_y_{i,j+\frac{1}{2}}^{(2k+1)}+E_y_{i,j+\frac{1}{2}}^{(2k)})
\end{equation}
\begin{equation}
\mu  \left[(H^{*}_z_{i+\frac{1}{2},j+\frac{1}{2}})^2 -(H_z_{i+\frac{1}{2},j+\frac{1}{2}}^{(2k)})^2 \right] = -\frac{\Delta t}{2}\varLambda_x {(E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}})} (H^{*}_z_{i+\frac{1}{2},j+\frac{1}{2}}+H_z_{i+\frac{1}{2},j+\frac{1}{2}}^{(2k)})
\end{equation} \\
Summing over all the terms in the above two equations, and add them together. Note that $E_y$ satisfies the boundary condition $(E,0) \times (\textbf{n},0) =0$ on (0,T] x $\partial \Omega$, where \textbf{n} is the outward normal vector on $\partial \Omega$, then by Lemma 9 in \cite{Chen2008} we have
\begin{equation}
\sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} (\epsilon(E{_y_{i,j+\frac{1}{2}}}^{(2k+1)})^2 + \mu(H_z_{i+\frac{1}{2},j+\frac{1}{2}}^*)^2 ) = \sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} (\epsilon(E{_y_{i,j+\frac{1}{2}}}^{(2k)})^2 + \mu(H_z_{i+\frac{1}{2},j+\frac{1}{2}}^{(2k)})^2 )
\end{equation} 
Similarly, from (2.9), we have,
\begin{equation}
\sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} (\epsilon(E{_x_{i+\frac{1}{2},j}}^{(2k+1)})^2 + \mu(H_z_{i+\frac{1}{2},j+\frac{1}{2}}^*)^2 ) = \sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} (\epsilon(E{_x_{i+\frac{1}{2},j}}^{(2k)})^2 + \mu(H_z_{i+\frac{1}{2},j+\frac{1}{2}}^{*})^2 )
\end{equation} 

Similarly, the same can be done for the even time step.\\
Equations (2.42) and (2.43) satisfy the condition in equation (2.38). \\
For (2.39), let us consider one time period (n = 2k), we use (2.11) and (2.13). Subtracting the (2k)th step from (2k+2)th step of the EO-4th-EC scheme, we get:
\begin{multline}
\frac{ E^{(2k+2)}_{x_{i+\frac{1}{2},j}}-E^{(2k)}_{x_{i+\frac{1}{2},j}} } {\varDelta t} - \frac{ E^{(2k+1)}_{x_{i+\frac{1}{2},j}}-E^{(2k-1)}_{x_{i+\frac{1}{2},j}} } {\varDelta t} 
\\=
  \frac{1}{2\epsilon} \varLambda_y \{ (H^{**}_{z_{i+\frac{1}{2},j}} -  H^{**-2}_{z_{i+\frac{1}{2},j}}) + (H^{(2k+1)}_{z_{i+\frac{1}{2},j}}-H^{(2k-1)}_{z_{i+\frac{1}{2},j}} )\}, 
\end{multline}
\begin{multline}
\frac{H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{**-2}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} -\frac{H^{2k+1}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k-1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} \\= 
\frac{1}{2\mu}\varLambda_y\{ (E^{(2k+2)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}-E^{(2k)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}})+( E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}-E^{(2k-1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}})\},
\end{multline}
\begin{multline}
\frac{ E^{(2k+2)}_{y_{i,j+\frac{1}{2}}}-E^{(2k)}_{y_{i,j+\frac{1}{2}}} } {\varDelta t} - \frac{ E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}-E^{(2k-1)}_{y_{i,j+\frac{1}{2}}} } {\varDelta t} 
\\=
-\frac{1}{2\epsilon} \varLambda_x \{ (H^{(2k+2)}_{z_{i,j+\frac{1}{2}}} -  H^{(2k)}_{z_{i,j+\frac{1}{2}}}) + (H^{**}_{z_{i,j+\frac{1}{2}}}-H^{**-2}_{z_{i,j+\frac{1}{2}}} )\}, 
\end{multline}
\begin{multline}
\frac{H^{(2k+2)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} -\frac{H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(**-2)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} \\= 
\frac{1}{2\mu}\varLambda_x\{ (E^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}-E^{(2k)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}})+( E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}-E^{(2k-1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}})\},
\end{multline}

Note that $\frac{ E^{(2k+2)}-E^{(2k)} } {2\varDelta t}$ also satisfies the boundary condition, then following the proof of  discrete energy conservation I and using Lemma 9 in paper \cite{Chen2008}, we can get the second discrete energy conservation. 
\end{proof}

Theorem 2.5.1 shows that the EO-4th-EC scheme satisfies the Poynting theorem in discrete sense and the energy conservation in discrete time variation holds at every two time steps.\\

\section{Convergence Analysis}
\begin{lemma}[Truncation Error]
%	{Lemma: Truncation Error}
\textit{Assume that the solutions are smooth enough, for example \textbf{E} and $H\in \textbf{C}^5$ $([0,T];[\textbf{C}^5(\Omega)]^2)$. Then the truncation errors are fourth order in space and second order in time.}
\begin{equation}
max{\mid \xi_x^{n+\frac{1}{2}} \mid,\mid \xi_y^{n+\frac{1}{2}}\mid, \mid \eta_z^{n+\frac{1}{2}}\mid} \leqslant C_2 (\epsilon, \mu)\left({\Delta t^2+ \Delta x^4+\Delta y^4 }\right), 
\end{equation}
Where $C_2(\epsilon ,\mu)$ is constant dependent of parameter $\epsilon$ and $\mu$.\\	
\end{lemma}
	


Now, we provide the convergence analysis for the scheme. Let $\xi_{x_{i+\frac{1}{2},j}}^{2k+1}, \xi_{y_{i,j+\frac{1}{2}}}^{2k+1}$ and  $\eta_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{2k+1}$ be the truncation error of at time $t^{2k+1}$.

Now, we will provide an analysis of convergence of the EO-4th-EC scheme.

\begin{theorem}[Convergence Analysis]
\textit{Suppose that the exact solutions of the Maxwell's equation $E_x$, $E_y$ and $H_z$ are smooth enough:}
\begin{equation}
\vec{E}\in C^3 ([0,T];[C^3(\Omega)^2]) and H_z \in C^3([0,T];C^3(\Omega)).
\end{equation}

Let $E_x^n$, $E_y^n$ and $H_z^n$ be the numerical solutions of the EO-4th-EC scheme. Then, for any fixed time interval $T>0$, there exists a constant $C_{\epsilon\mu}$ independent of $\Delta t$, $\Delta x$ and $\Delta y$ such that $n \geqslant0$,
\begin{multline}
\max\limits_{0 \leq n \leq N} { \|\epsilon^{\frac{1}{2}} [E(t^n)-E^n]\|^2_E +  \|\mu^{\frac{1}{2}} [H_z(t^n)-H^n_z\|^2_H ]}  \\ 
\leq e^T ( \|\epsilon^{\frac{1}{2}} [E(t^0)-E^0]\|^2_E +  \|\mu^{\frac{1}{2}} [H_z(t^0)-H^0_z\|^2_H]) +C_{\mu\epsilon}e^T (\Delta t^2 + \Delta x^4 +\Delta y^4)^2
\end{multline}
where n=2k and 2k+1.
\end{theorem}

\begin{proof}
	 Firstly, let us define the error functions on the staggered grid:
\begin{equation}
\varepsilon_{\omega_{\alpha,\beta}}^n = E_\omega(x_\alpha ,y_\beta,t^n) - E_{\omega_{\alpha,\beta}}^n, H _{z_{\alpha,\beta}}^n = H_z(x_\alpha ,y_\beta ,t^n)-H_{z_{\alpha,\beta}}^n,
\end{equation}

where $E_\omega(x_\alpha ,y_\beta ,t^n) (\omega=x,y)$ and $H_z(x_\alpha ,y_\beta ,t^n)$ denote the values of exact solutions $E_\omega$ and $H_z$ on the point $(x_\alpha ,y_\beta ,t^n)$. Then, we focus on the following error equations:
\begin{equation}
\frac{\varepsilon^{(2k+1)}_{y_{i,j+\frac{1}{2}}}-\varepsilon^{(2k)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t} =- \frac{1}{2\epsilon}\varLambda_x\{H^*_{z_{i,j+\frac{1}{2}}} + H^{(2k)}_{z_{i,j+\frac{1}{2}}}\} + e_1_{i,j+\frac{1}{2}}, 
\end{equation}
\begin{equation}
\frac{H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t}
= - \frac{1}{2\mu}\varLambda_x\{\varepsilon^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + \varepsilon^{(2k)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\}  + e_2_{i+\frac{1}{2},j+\frac{1}{2}},
\end{equation}
\begin{equation}
\frac{\varepsilon^{(2k+1)}_{x_{i+\frac{1}{2},j}}-\varepsilon^{(2k)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} =  \frac{1}{2\epsilon}\varLambda_y\{H^{(2k+1)}_{z_{i+\frac{1}{2},j}} + H^{*}_{z_{i+\frac{1}{2},j}}\}  + e_3_{i+\frac{1}{2},j}, 
\end{equation}
\begin{equation}
\frac{H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu} \varLambda_y \{\varepsilon^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + \varepsilon^{(2k)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\}  + e_4_{i+\frac{1}{2},j+\frac{1}{2}},
\end{equation}
\begin{equation}
\frac{\varepsilon^{(2k+2)}_{x_{i+\frac{1}{2},j}}-\varepsilon^{(2k+1)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} =  \frac{1}{2\epsilon}\varLambda_y\{H^{**}_{z_{i+\frac{1}{2},j}} + H^{(2k+1)}_{z_{i+\frac{1}{2},j}}\}  + e_5_{i,j+\frac{1}{2}}, 
\end{equation}
\begin{equation}
\frac{H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu}\varLambda_y\{\varepsilon^{(2k+2)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + \varepsilon^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\}  + e_6_{i+\frac{1}{2},j+\frac{1}{2}},
\end{equation}
\begin{equation}
\frac{\varepsilon^{(2k+2)}_{y_{i,j+\frac{1}{2}}}-\varepsilon^{(2k+1)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t} =  -\frac{1}{2\epsilon}\varLambda_x\{H^{(2k+2)}_{z_{i,j+\frac{1}{2}}} + H^{**}_{z_{i,j+\frac{1}{2}}}\} + e_7_{i,j+\frac{1}{2}}, 
\end{equation}
\begin{equation}
\frac{H^{(2k+2)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} =-\frac{1}{2\mu}\varLambda_x\{\varepsilon^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + \varepsilon^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\} + e_8_{i+\frac{1}{2},j+\frac{1}{2}},
\end{equation} 

where $e_1_{i,j+\frac{1}{2}}$, $ e_2_{i+\frac{1}{2},j}$,  $e_3_{i+\frac{1}{2},j+\frac{1}{2}}$, etc.. are truncation errors.\\ 

Without loss of generality, we regard $\varepsilon^{(2k+1)}_{x_{i+\frac{1}{2},j}}$, $\varepsilon^{(2k+1)}_{y_{i,j+\frac{1}{2}}}$, $H^{*}_{z_{i+\frac{1}{2},j}}$ and $H^{**}_{z_{i+\frac{1}{2},j}}$ as the temperate variable from (2k+1)-th step to the (2k+2)-th step. Thus we get:\\$ $\\
$e_1_{i,j+\frac{1}{2}} = \frac{1}{2} \xi_{y_{i,j+\frac{1}{2}}}^{(2k+1)}  $, $e_2_{i+\frac{1}{2},j+\frac{1}{2}} = 0$, $e_3_{i+\frac{1}{2},j} = \frac{1}{2} \xi_{x_{i+\frac{1}{2},j}}^{(2k+1)}$, $e_4_{i+\frac{1}{2},j+\frac{1}{2}} = \frac{1}{2} \eta_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{(2k+1)}$, \\$ $\\
$e_5_{i+\frac{1}{2},j} = \frac{1}{2} \xi_{x_{i+\frac{1}{2},j}}^{(2k+1)}  $, $e_6_{i+\frac{1}{2},j+\frac{1}{2}} = 0$, 
$e_7_{i,j+\frac{1}{2}} = \frac{1}{2} \xi_{y_{i,j+\frac{1}{2}}}^{(2k+1)}$, $e_8_{i+\frac{1}{2},j+\frac{1}{2}} = \frac{1}{2} \eta_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{(2k+1)}$, \\ $ $\\

Note that the truncation errors $\xi_{y_{i,j+\frac{1}{2}}}^{(2k+1)}$, $\xi_{x_{i+\frac{1}{2},j}}^{(2k+1)}$ and $\eta_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{(2k+1)}$ satisfy  (2.46). From (2.52)-(2.53) and using Lemma 9 in \cite{Chen2008}, we get:
\begin{multline}
\sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} (\epsilon(\varepsilon_{y_{i,j+\frac{1}{2}}}^{(2k+1)})^2 + \mu(H_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^*)^2 - \epsilon(\varepsilon_{y_{i,j+\frac{1}{2}}}^{(2k)})^2 - \mu(H_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{(2k)})^2) \\
= \frac{1}{2} \Delta t
\sum\limits_{i=0}^{I-1} \sum\limits_{j=0}^{J-1} (\epsilon\xi_{y_{i,j+\frac{1}{2}}}^{(2k+1)}(\varepsilon_{y_{i,j+\frac{1}{2}}}^{(2k+1)}+\varepsilon_{y_{i,j+\frac{1}{2}}}^{(2k)}) +
 \mu\eta_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{(2k+1)}(H_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^*+H_{z_{i+\frac{1}{2},j+\frac{1}{2}}}^{(2k)}))
\end{multline}

By the Schwartz inequality, we have
\begin{multline}
(1-\frac{\Delta t}{4})  { \epsilon\|[\varepsilon_y^{2k+1} \|^2_{E_y} +
 	\mu \| H^*_z\|^2_H }  \\ 
\leq C_1 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
(1+\frac{\Delta t}{4})  ( \epsilon\|\varepsilon_y^{2k} \|^2_{E_y} +
	\mu \| H^{(2k)}_z\|^2_H ).
\end{multline} 

Similarly it follows from (2.54) to (2.59) that
\begin{multline}
(1-\frac{\Delta t}{4})  { \epsilon\|[\varepsilon_x^{(2k+1)} \|^2_{E_x} + \mu\| H^{(2k+1)}_z\|^2_H }  \\ 
\leq C_1\Delta t (\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
(1+\frac{\Delta t}{4})  ( \epsilon\|\varepsilon_x^{(2k)} \|^2_{E_x} +
\mu \| H^{*}_z\|^2_H ),
\end{multline} 
\begin{multline}
(1-\frac{\Delta t}{4})  { \epsilon\|[\varepsilon_x^{(2k+2)} \|^2_{E_x} + \mu\| H^{**}_z\|^2_H }  \\ 
\leq C_1 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
(1+\frac{\Delta t}{4})  ( \epsilon\|\varepsilon_x^{2k+1} \|^2_{E_x} +
\mu \| H^{(2k+1)}_z\|^2_H ),
\end{multline} 
and
\begin{multline}
(1-\frac{\Delta t}{4})  { \epsilon\|[\varepsilon_y^{2k+2} \|^2_{E_y} +
	\mu \| H^{(2k+2)}_z\|^2_H }  \\ 
\leq C_1 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
(1+\frac{\Delta t}{4})  ( \epsilon\|\varepsilon_y^{(2k+1)} \|^2_{E_y} +
\mu \| H^{**}_z\|^2_H ).
\end{multline} 
Dividing both sides of (2.61) and (2.63) with $(1-\Delta t/4)$, both sides of (2.62) and (2.64) with $(1+\Delta t/4)$, and adding (2.61) to (2.62) and (2.63) to (2.64) respectively, we can get
\begin{multline}
\epsilon\|\varepsilon_y^{(2k+1)} \|^2_{E_y} + \frac{(1-\frac{\Delta t}{4})}{(1+\frac{\Delta t}{4})}  ( \epsilon\|[\varepsilon_x^{2k+1} \|^2_{E_x} + \mu \| H^{(2k+1)}_z\|^2_H )  \\ 
\leq C_2 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
\epsilon\|\varepsilon_x^{(2k)} \|^2_{E_x} +
\frac{(1+\frac{\Delta t}{4})}{(1-\frac{\Delta t}{4})}   ( \epsilon\|\varepsilon_y^{(2k)} \|^2_{E_y} +
\mu \| H^{(2k)}_z\|^2_H ).
\end{multline} 
\begin{multline}
\epsilon\|\varepsilon_x^{(2k+2)} \|^2_{E_x} + \frac{(1-\frac{\Delta t}{4})}{(1+\frac{\Delta t}{4})}  ( \epsilon\|[\varepsilon_y^{2k+2} \|^2_{E_y} + \mu \| H^{(2k+2)}_z\|^2_H )  \\ 
\leq C_3 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
\epsilon\|\varepsilon_y^{(2k+1)} \|^2_{E_y} +
\frac{(1+\frac{\Delta t}{4})}{(1-\frac{\Delta t}{4})}   ( \epsilon\|\varepsilon_x^{(2k+1)} \|^2_{E_x} +
\mu \| H^{(2k+1)}_z\|^2_H ).
\end{multline} 
Multiplying both sides of (2.66) with $(1-\frac{\Delta t}{4})^2(1+\frac{\Delta t}{4})^{-2}$, we get:
\begin{multline}
\frac{(1-\frac{\Delta t}{4})^2}{(1+\frac{\Delta t}{4})^{2}} \epsilon\|\varepsilon_x^{(2k+2)} \|^2_{E_x} + \frac{(1-\frac{\Delta t}{4})^3}{(1+\frac{\Delta t}{4})^3}  ( \epsilon\|[\varepsilon_y^{2k+2} \|^2_{E_y} + \mu \| H^{(2k+2)}_z\|^2_H )  \\ 
\leq C_4 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
\frac{(1-\frac{\Delta t}{4})^2}{(1+\frac{\Delta t}{4})^{2}}\epsilon\|\varepsilon_y^{(2k+1)} \|^2_{E_y} +
\frac{(1-\frac{\Delta t}{4})}{(1+\frac{\Delta t}{4})}   ( \epsilon\|\varepsilon_x^{(2k+1)} \|^2_{E_x} +
\mu \| H^{(2k+1)}_z\|^2_H )\\
\leq C_4 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
\epsilon\|\varepsilon_y^{(2k+1)} \|^2_{E_y} +
\frac{(1-\frac{\Delta t}{4})}{(1+\frac{\Delta t}{4})}   ( \epsilon\|\varepsilon_x^{(2k+1)} \|^2_{E_x} +
\mu \| H^{(2k+1)}_z\|^2_H )
\end{multline}

Adding (2.65) and (2.67) gives us
\begin{multline}
\frac{(1-\frac{\Delta t}{4})^2}{(1+\frac{\Delta t}{4})^{2}} \epsilon\|\varepsilon_x^{(2k+2)} \|^2_{E_x} + \frac{(1-\frac{\Delta t}{4})^3}{(1+\frac{\Delta t}{4})^3}  ( \epsilon\|[\varepsilon_y^{2k+2} \|^2_{E_y} + \mu \| H^{(2k+2)}_z\|^2_H )  \\ 
\leq C_5 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
\epsilon\|\varepsilon_x^{(2k)} \|^2_{E_x} +
\frac{(1+\frac{\Delta t}{4})}{(1-\frac{\Delta t}{4})}   ( \epsilon\|\varepsilon_y^{(2k)} \|^2_{E_y} +
\mu \| H^{(2k+1)}_z\|^2_H )
\end{multline}

Note that $\frac{(1+\frac{\Delta t}{4})}{(1-\frac{\Delta t}{4})} >1$. Therefore,
\begin{multline}
 \epsilon\|[\varepsilon_x^{2k+2} \|^2_{E_x} + \epsilon\|[\varepsilon_y^{2k+2} \|^2_{E_y} + \mu \| H^{(2k+2)}_z\|^2_H \\
 \leq C_6 \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
 \frac{(1+\frac{\Delta t}{4})^4}{(1-\frac{\Delta t}{4})^4} (\epsilon\|\varepsilon_x^{(2k)} \|^2_{E_x} + \epsilon\|\varepsilon_y^{(2k)} \|^2_{E_y} +
 \mu \| H^{(2k)}_z\|^2_H )
\end{multline}

Since $\Delta t=T/N$ and the total number of time steps N is even, i.e., N=2m, it can be verified that
$$ \left( \frac{(1+\frac{\Delta t}{4})^4}{(1-\frac{\Delta t}{4})^4} \right)^m \leq e^T. $$

Computing the inequality (2.69) from the (2k+2)-th step to the initial time step recursively by every two neighbouring steps, we get
\begin{multline}
\epsilon\|[\varepsilon_x^{2k+2} \|^2_{E_x} + \epsilon\|[\varepsilon_y^{2k+2} \|^2_{E_y} + \mu \| H^{(2k+2)}_z\|^2_H \\
\leq C_7 e^T \Delta t(\Delta t^2 + \Delta x^4 +\Delta y^4)^2 + 
e^T (\epsilon\|\varepsilon_x^{(0)} \|^2_{E_x} + \epsilon\|\varepsilon_y^{(0)} \|^2_{E_y} +
\mu \| H^{(0)}_z\|^2_H )
\end{multline}

Note that from the inequality (2.65), the error at the odd step can be controlled by that at the even step. Then the theorem is proved.
\end{proof}


From the Maxwell's equations, the electromagnetic waves must be divergence free since we assume that the electric charge density is zero: $\rho=0$, i.e.,
\begin{equation}
div(\epsilon\vec{E})=0.
\end{equation}

The following result demonstrates that the divergence free constraint is satisfied approximately.
\begin{theorem}[Convergence of the discrete divergence free]
\textit{Consider the EO-4th-EC scheme and assume that the total number of total time steps N is even, N=2m. Then we have:}
\begin{multline}
\epsilon(\Lambda_x E_{x_{i,j}}^{(2k+1)}+\Lambda_y E_{y_{i,j}}^{(2k+1)}+\Lambda_x E_{x_{i,j}}^{(2k+2)}+\Lambda_y E_{y_{i,j}}^{(2k+2)})+ \frac{\Delta t^2}{2}\Lambda_x \Lambda_y \delta_t H_{z_{i,j}}^{2k+1+\frac{1}{2}} \\
= \epsilon(\Lambda_x E_{x_{i,j}}^{(0)}+\Lambda_y E_{y_{i,j}}^{(0)}+\Lambda_x E_{x_{i,j}}^{(1)}+\Lambda_y E_{y_{i,j}}^{(1)})+ \frac{\Delta t^2}{2}\Lambda_x \Lambda_y \delta_t H_{z_{i,j}}^{\frac{1}{2}}.
\end{multline}
\end{theorem}
\begin{proof}
From equations (2.17)-(2.22) and noting that $\Lambda_x$ (or$\Lambda_y$) and $\delta_t$ are changeable, we get
\begin{equation} 
\delta_t \Lambda_x E^{(2k+\frac{1}{2})}_{x_{i,j}}  =  {\frac{1}{2\epsilon}} {\Lambda_x \varLambda_y}\left\{H^{(2k+1)}_{z_{i,j}} + H^{(2k)}_{z_{i,j}}\right\} -\frac{\Delta t}{4\mu\epsilon}{\varLambda_x}{\Lambda_x}{\varLambda_y}\left\{E^{(2k+1)}_{y_{i,j}}+E^{(2k)}_{y_{i,j}}\right\},
\end{equation}
\begin{equation}
\delta_t \Lambda_y E^{(2k+\frac{1}{2})}_{y_{i,j}}  = -\frac{1}{2\epsilon} {\varLambda_x}{\Lambda_y}\left\{H^{(2k+1)}_{z_{i,j}} + H^{(2k)}_{z_{i,j}}\right\} 
+\frac{\Delta t}{4\mu\epsilon}\Lambda_y\varLambda_x\varLambda_y \left\{E^{(2k+1)}_{x_{i,j}}+E^{(2k)}_{x_{i,j}}\right\},
\end{equation}
\begin{equation} 
\delta_t \Lambda_x E^{(2k+1+\frac{1}{2})}_{x_{i,j}}  =  {\frac{1}{2\epsilon}} {\Lambda_x \varLambda_y}\left\{H^{(2k+2)}_{z_{i,j}} + H^{(2k+1)}_{z_{i,j}}\right\} +\frac{\Delta t}{4\mu\epsilon}{\varLambda_x}{\Lambda_x}{\varLambda_y}\left\{E^{(2k+2)}_{y_{i,j}}+E^{(2k+1)}_{y_{i,j}}\right\},
\end{equation}
\begin{equation}
\delta_t \Lambda_y E^{(2k+1+\frac{1}{2})}_{y_{i,j}}  = -\frac{1}{2\epsilon} {\varLambda_x}{\Lambda_y}\left\{H^{(2k+2)}_{z_{i,j}} + H^{(2k+1)}_{z_{i,j}}\right\} 
-\frac{\Delta t}{4\mu\epsilon}\Lambda_y\varLambda_x\varLambda_y \left\{E^{(2k+2)}_{x_{i,j}}+E^{(2k+1)}_{x_{i,j}}\right\},
\end{equation}
Summing up (2.73) to (2.76), we get:
\begin{multline}
\delta_t \Lambda_x E^{(2k+\frac{1}{2})}_{x_{i,j}} + \delta_t \Lambda_y E^{(2k+\frac{1}{2})}_{y_{i,j}} + \delta_t \Lambda_x E^{(2k+1+\frac{1}{2})}_{x_{i,j}} + \delta_t \Lambda_y E^{(2k+1+\frac{1}{2})}_{y_{i,j}} \\
= \frac{\Delta t}{2\epsilon}\varLambda_x\varLambda_y \left\{ \frac{1}{2\mu} \Lambda_y(E^{(2k+1)}_{x_{i,j}}+E^{(2k)}_{x_{i,j}}) -\frac{1}{2\mu} \Lambda_x(E^{(2k+1)}_{y_{i,j}}+E^{(2k)}_{y_{i,j}})\\
-\frac{1}{2\mu} \Lambda_y(E^{(2k+2)}_{x_{i,j}}+E^{(2k+1)}_{x_{i,j}})
+\frac{1}{2\mu} \Lambda_x(E^{(2k+2)}_{y_{i,j}}+E^{(2k+1)}_{x_{i,j}})\right\} \\
= \frac{\Delta t}{2\epsilon}\varLambda_x\varLambda_y \left( \delta_t H_{z_{i,j}}^{(2k+\frac{1}{2})} - \delta_t H_{z_{i,j}}^{(2k+1+\frac{1}{2})} \right) = -\frac{\Delta t^2}{2\epsilon}\Lambda_x \Lambda_y \delta_t \delta_t H_{z_{i,j}}^{(2k+1)}.
\end{multline}
Consequently, we have
\begin{multline}
\Lambda_x E^{(2k+1)}_{x_{i,j}} + \Lambda_y E^{(2k+1)}_{y_{i,j}} +\Lambda_x E^{(2k+2)}_{x_{i,j}} + \Lambda_y E^{(2k+2)}_{y_{i,j}} + \frac{\Delta t^2}{2\epsilon}\Lambda_x \Lambda_y \delta_t  H_{z_{i,j}}^{(2k+1+\frac{1}{2})}\\
=\Lambda_x E^{(0)}_{x_{i,j}} + \Lambda_y E^{(0)}_{y_{i,j}} +\Lambda_x E^{(1)}_{x_{i,j}} + \Lambda_y E^{(1)}_{y_{i,j}} + \frac{\Delta t^2}{2\epsilon}\Lambda_x \Lambda_y \delta_t  H_{z_{i,j}}^{(\frac{1}{2})}
\end{multline}
This completes the proof of the theorem. 
\end{proof}


\section{Numerical Results}
\label{sec:num}
\subsection{Phase Error}
We can now present the numerical dispersion errors of the proposed scheme and other schemes with different wave propagation angles, different grid sizes and CFL numbers. \\

Let $\xi=e^{(i\omega 2 \Delta t)}$, where $\omega$ is a complex number: $\omega=\omega_R + \omega_I$ are the real and imaginary parts respectively. Then,
\begin{equation}
\xi = e^{-\omega_I 2\Delta t}(\cos(\omega_R 2\Delta t)+i\sin(\omega_R 2\Delta t)).
\end{equation}
Let $Im(\xi)$ and $Re(\xi)$ be the imaginary and real parts of the growth factor, $\xi$. The numerical phase velocity, $v_p$ ,is normalized to the speed of the light, c, also called phase velocity error, as:
\begin{equation}
\frac{v_p}{c}= \frac{N_\lambda}{2\pi S}\arctan{\left(\frac{|Im(\xi)|}{|Re(\xi)|}\right)}
\end{equation}

where
wavenumber \textbf{k} = $(k_x,k_y)$ : $k_x = k \cos{\phi}$, $k_y = k \sin{\phi}$,
$\omega = ck$, $c=\frac{1}{\sqrt{\mu\epsilon}}$,
$\lambda $= wavelength,
$\Delta x$ = $\Delta y$ = h  be the spatial step size,
$N_{\lambda} = \frac{\lambda}{h} $ = number of points per wavelength (NPPW),
$S=\frac{c \Delta t}{h}$ = CFL (Courant) number.\\


\begin{figure}
\includegraphics[width=0.8\textwidth]{Figure13}\\
\caption{ Numerical dispersion against the wave propagation angle phi from 0 to 180 degrees with $N_{\lambda}$=40 and S=3.5.}
\end{figure}

\begin{figure}
\includegraphics[width=0.8\textwidth]{Figure2}\\
\caption{ Numerical dispersion against number of points per wavelength, $N_{\lambda}$ with S=3.5, $\phi$=65 degrees.}
\end{figure}

\begin{figure}
\includegraphics[width=0.8\textwidth]{Figure1}\\
\caption{ Numerical dispersion against CFL number  with $N_{\lambda}$=40 and $\phi$=65 degrees.} 
\end{figure}
For our even-odd cycled energy conserved scheme, equation (2.80) is used for which the real and imaginary parts of the stability factors are used. The analytical phase error of the Maxwell's equations is 1. Figure 1 shows the numerical phase error against wave propagation angle, $\phi$, as it varies from $0^o$ to $180^o$ using parameters $N_{\lambda}$=40 and S=3.5. Figure 2 compares the numerical phase error against number of points per wavelength as they vary from 10 to 100 with parameters $N_{\lambda}$ with S=3.5, $\phi =65^o$. Figure 3 compares the phase error against Courant number,CFL, with parameters $N_{\lambda}$=40 and $\phi = 65^o$, where CFL number varies from 1 to 5. All three figures compare the proposed scheme with Crank Nicolson(CN), EC-S-FDTDII [1], ADI-FDTD [7,13], EC-S-FDTD(1,4) and EC-S-FDTD(2,4) [3] schemes. As shown in Figures 1, 2 and 3, the CN scheme has the worst phase velocity error comparing to other schemes. ADI-FDTD and EC-S-FDTDII schemes are although better than the CN scheme, they are not as great as the EC-S-FDTD(2,4) or EC-S-FDTD(1,4) scheme. Comparing our proposed scheme to all the other schemes, our scheme has the ratio $v_p/c$ closest to 1. Also, Figure 3 shows that as CFL number increases, the phase error increases, whereas in Figure 2, as the number of points per wavelength increases, the phase velocity error decreases. Therefore, the numerical dispersion error of our scheme, EO-4th-EC scheme has the smallest error comparing to the other schemes.


\singlespacing

\subsection{Accuracy, Energy Conservation and Divergence Free Approximation}
\doublespacing
Now, we will focus on energy conservation, accuracy and convergence of divergence free by comparing the EC-S-FDTD(1,4) (named EC(1,4)), EC-S-FDTD(2,4) (named EC(2,4)) [3], ADI-FDTD [7, 13] schemes to our EO-4th-EC scheme. We consider the Maxwell's equations (1)-(3) in a lossless medium with the domain surrounded by a perfect conductor. The exact solution to Eqs. (1) - (3) is
\begin{equation}
E_x=\frac{k_y}{\epsilon \sqrt{\mu}\omega} \cos{(\omega \pi t)} \cos{(k_x \pi x)} \sin{(k_y\pi y)}
\end{equation}
\begin{equation}
E_y=-\frac{k_x}{\epsilon \sqrt{\mu}\omega} \cos{(\omega \pi t)} \sin{(k_x \pi x)} \cos{(k_y\pi y)}
\end{equation}
\begin{equation}
H_z=-\frac{1}{\sqrt{\mu}} \sin{(\omega \pi t)} \cos{(k_x \pi x)} \cos{(k_y\pi y)}
\end{equation}

The exact energy of electromagnetic fields can be computed directly as $$\text{EnergyI} = \left( \int_{\Omega}(\epsilon|\textbf{E}(x,y,t)|^2+\mu|H_z (x,y,t)|^2)dxdy \right)^\frac{1}{2} =\frac{1}{2}.$$
Define the relative errors of energy by:
\begin{equation}
\text{REE-I}=\max_{\limits{0 \leq n \leq N}} \frac{|(\|\epsilon^{\frac{1}{2}} \textbf{E}^n\|^2_E +(\|\mu^{\frac{1}{2}} H_z^n\|^2_H  )^{\frac{1}{2}}-EnergyI |}{EnergyI}
\end{equation}
and
\begin{equation}
\text{REE-II}=\max_{\limits{0 \leq n \leq N-1}} \frac{|(\|\epsilon^{\frac{1}{2}} \delta_t\textbf{E}^{n+\frac{1}{2}}\|^2_E +(\|\mu^{\frac{1}{2}} \delta_t H_z^{n+\frac{1}{2}}\|^2_H  )^{\frac{1}{2}}
	- (\|\epsilon^{\frac{1}{2}} \delta_t\textbf{E}^{\frac{1}{2}}\|^2_E +(\|\mu^{\frac{1}{2}} \delta_t H_z^{\frac{1}{2}}\|^2_H  )^{\frac{1}{2}}
	|}{(\|\epsilon^{\frac{1}{2}} \delta_t\textbf{E}^{\frac{1}{2}}\|^2_E +(\|\mu^{\frac{1}{2}} \delta_t H_z^{\frac{1}{2}}\|^2_H  )^{\frac{1}{2}}}
\end{equation}

Let ErrorI and ErrorII be defined by:
\begin{equation}
\text{ErrorI}=\max_{\limits{0 \leq n \leq N}} \frac{(\|\epsilon^{\frac{1}{2}} [\textbf{E}(t^n)-E^n]\|^2_E +\|\mu^{\frac{1}{2}} [H_z(t^n) -H_z^n]\|^2_H  )^{\frac{1}{2}}}{EnergyI}
\end{equation}
and
\begin{equation}
\text{ErrorII}=\max_{\limits{0 \leq n \leq N-1}} \frac{|(\|\epsilon^{\frac{1}{2}} [\delta_t\textbf{E}(t^n) -\delta_t E^n ]\|^2_E +(\|\mu^{\frac{1}{2}} [\delta_t H_z(t^{n})-\delta_t H_z^n   ]\|^2_H  )^{\frac{1}{2}}} {EnergyII}
\end{equation}
\begin{table}	
\caption{Relative Errors of Energy I and Energy II of different schemes, Parameters: T = 1, $\Delta$t = $\Delta$x = $\Delta$y = T/100 and $k_x = k_y=k$} \\
\small
\begin{tabular}{|l|l|l|l|l|l|l|l|l|} 
\hline
 $k$ & \multicolumn{2}{c|}{EC(1,4)} & \multicolumn{2}{c|}{EC(2,4)} & \multicolumn{2}{c|}{EO-4th-EC} & \multicolumn{2}{c|}{ADI-FDTD}   \\
\hline
 & REE-I & REE-II & REE-I & REE-II & REE-I & REE-II & REE-I & REE-II  \\
\hline
1& 2.021e-14 & 1.760e-14	& 1.399e-14 &	 1.260e-14	 &1.976e-14	 &1.780e-14 & 6.165e-05 & 2.876e-04
\\
\hline
3&	1.976e-14&	1.856e-14	&1.399e-14	 &1.268e-14	 &1.976e-14	 &1.807e-14 & 5.532e-04 & 0.0026
\\
\hline
7	&1.932e-14	&1.937e-14	&1.421e-14&	 1.416e-14	 &1.932e-14	 &1.942e-14 & 0.003 & 0.0139
\\
\hline
\end{tabular}
\end{table}
\normalsize
\begin{table*}
\centering
\caption{Relative Error I in time step of numerical solutions of difference schemes, Parameters: T=1, $k_x = k_y$ =1, $\mu=\epsilon$=1, $\Delta t= \Delta x= \Delta y = 1/N$}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
$N$ &  \multicolumn{2}{c|}{EC(1,4)} &  \multicolumn{2}{c|}{EC(2,4)} &  \multicolumn{2}{c|}{EO-4th-EC} & \multicolumn{2}{c|}{ADI-FDTD} \\
\hline
 & Error I & Ratio &Error I & Ratio & Error I & Ratio & Error I & Ratio\\
\hline
25&	0.0434	&--	&0.0051&	--&	0.0126&	-- & 0.0108 & -- \\
\hline
50&	0.0215&	1.0134	&0.0013	&1.9854&	0.0032	 &1.9773 & 0.0027 & 2.0000 \\
\hline
100	&0.0107&	1.0067	&3.2209e-04	&1.9999	 &7.9476e-04	 &2.0095  &6.7599e-04 & 1.9979\\
\hline
200	&0.0054	&0.9866&	8.0527e-05&	1.9999&	 1.9870e-04	 &1.9999 & 1.6902e-04 & 1.9998 \\
\hline
\end{tabular}
\end{table*}
\begin{table*}
\centering
\caption{Relative Error II in time step of numerical solutions of difference schemes, Parameters: T=1, $k_x = k_y$ =1, $\mu=\epsilon$=1, $\Delta t= \Delta x= \Delta y = 1/N$}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
$N$ & \multicolumn{2}{c|}{EC(1,4)} & \multicolumn{2}{c|}{EC(2,4)}  & \multicolumn{2}{c|}{EO-4th-EC} & \multicolumn{2}{c|}{ADI-FDTD} \\
\hline
 & Error II & Ratio &Error II & Ratio & Error II & Ratio & Error II & Ratio\\
\hline
25&	0.0173	&--	&0.0051&	--&	0.0114&	-- & 0.0101 & --\\
\hline
50&	0.0071&	1.2849&0.0013&1.9720&	0.0031	 &1.8787 &0.0026  &1.9579\\
\hline
100	&0.0032&1.1497&3.2745e-04&1.9892&7.7781e-04	 &1.9948 & 6.4756e-04 & 2.0054 \\
\hline
200	&0.0015	&1.0931&8.2070e-05&	1.9963&	1.9550e-04	 &1.9922 & 1.6246e-04 & 1.9949\\
\hline
\end{tabular}
\end{table*}
\begin{table*}
\centering
\caption{Relative Error I in spatial step of numerical solutions of difference schemes, Parameters: T=1, $k_x = k_y$ =1, $\mu=\epsilon$=1, $\Delta t=\frac{1}{N^4}$ for EC(1,4) and $\Delta t= \frac{1}{N^2}$ for EC(2,4), ADI-FDTD and EO-4th-EC Scheme}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
$N$ & \multicolumn{2}{c|}{EC(1,4)} & \multicolumn{2}{c|}{EC(2,4)} &  \multicolumn{2}{c|}{EO-4th-EC} & \multicolumn{2}{c|}{ADI-FDTD} \\
\hline
 & Error I & Ratio & Error I & Ratio & Error I & Ratio & Error I & Ratio\\
\hline
10& 2.2837e-04  &  --    &	5.2279e-04	& --      & 9.8328e-04 & --     & 0.0187   & --       \\
\hline
15 & 4.5219e-05 & 3.9941  & 1.0341e-04	& 3.9966  & 1.9390e-04 & 4.0042 & 0.0082 & 2.0332  \\
\hline
 20& 1.4320e-05 & 3.9970  & 3.2735e-05	&3.9984	 & 6.1518e-05 & 3.9906 & 0.0046 & 2.0094  \\
\hline
 25& 5.8679e-06 & 3.9981  & 1.3411e-05	&3.9991	 & 2.5178e-05 & 4.0035 & 0.0029 & 2.0675  \\
\hline
 30& 2.8305e-06 & 3.9987  & 6.4682e-06	&3.9994	 & 1.2154e-05 & 3.9947 & 0.0020 & 2.0380  \\
\hline
\end{tabular}
\end{table*}
\begin{table*}
\centering
\caption{Relative Error II in spatial step of numerical solutions of difference schemes, Parameters: T=1, $k_x = k_y$ =1, $\mu=\epsilon$=1, $\Delta t=\frac{1}{N^4}$ for EC(1,4) and $\Delta t= \frac{1}{N^2}$ for EC(2,4), ADI-FDTD and EO-4th-EC Scheme}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
$N$ & \multicolumn{2}{c|}{EC(1,4)} &  \multicolumn{2}{c|}{EC(2,4)} & \multicolumn{2}{c|}{EO-4th-EC} & \multicolumn{2}{c|}{ADI-FDTD} \\
\hline
 & Error II & Ratio & Error II & Ratio & Error II & Ratio & Error II & Ratio\\
\hline
10& 2.0886e-04  &  --    &	5.3257e-04	& --      & 9.7475e-04 & --     & 0.0190   & --       \\
\hline
15 & 4.1383e-05 & 3.9941  & 1.0564e-04	& 3.9897  & 1.9298e-04 & 3.9944 & 0.0084 & 2.0130  \\
\hline
 20& 1.3108e-05 & 3.9970  & 3.3472e-05	&3.9951	 & 6.1464-05 & 3.9771 & 0.0047 & 2.0184  \\
\hline
 25& 5.3719e-06 & 3.9981  & 1.3719e-05	&3.9971	 & 2.5162e-05 & 4.0024 & 0.0030 & 2.0119  \\
\hline
 30& 2.5914e-06 & 3.9987  & 6.6182e-06	&3.9982	 & 1.2160e-05 & 3.9895 & 0.0021 & 1.9563  \\
\hline
\end{tabular}
\end{table*}
\begin{table*}
\centering
\caption{Error of divergence free of divergence 1 in time step of numerical solutions of different schemes, Parameters: T=1, $k_x = k_y$ =1, $\mu=\epsilon$=1, $\Delta t =\Delta x = \Delta y = \frac{1}{N}$}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
$N$ & \multicolumn{2}{c|}{EC(1,4)}  & \multicolumn{2}{c|}{EC(2,4)} & \multicolumn{2}{c|}{EO-4th-EC} & \multicolumn{2}{c|}{ADI-FDTD} \\
\hline
 & Div 1 & Ratio & Div 1 & Ratio & Div 1 & Ratio  & Div 1 & Ratio\\
\hline
25&	0.1964	&--	&0.0044&	--&	0.0349&	-- & 0.0174 & --\\
\hline
50&	0.0986&	0.9941&0.0011&2.0000&	0.0088	 &1.9877 & 0.0044 & 1.9835 \\
\hline
100	&0.0493&0.9493&2.7404e-04&1.9999&0.0022	 &2.0000 & 0.0019 &2.0711\\
\hline
200	&0.0247	&0.9971&6.8513e-05&	1.9999&	5.4810e-04	 &2.0050 &0.0011 & 1.8998 \\
\hline
\end{tabular}
\end{table*}
\begin{table*}
\centering
\caption{Error of divergence free of divergence 2 in time step of numerical solutions of different scheme, Parameters: T=1, $k_x = k_y$ =1, $\mu=\epsilon$=1, $\Delta t =\Delta x = \Delta y = \frac{1}{N}$}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|l|}
\hline
$N$ & \multicolumn{2}{c|}{EC(1,4)} & \multicolumn{2}{c|}{EC(2,4)} & \multicolumn{2}{c|}{EO-4th-EC} & \multicolumn{2}{c|}{ADI-FDTD} \\
\hline
 & Div 2 & Ratio & Div 2 & Ratio & Div 2 & Ratio & Div 2 & Ratio\\
\hline
25&	0.0982	&--	&0.0022&	--&	0.0175&	-- & 0.0087 & --\\
\hline
50&	0.0493&	0.9941&5.4791e-04&2.0055&	0.0044	 &1.9918 & 0.0022 & 1.9835 \\
\hline
100	&0.0247&0.9971&1.3702e-04&1.9996&0.0011	 &2.0000 & 9.7365e-04 & 2.0104 \\
\hline
200	&0.0123	&1.0059&3.4256e-05&	2.0000&	2.7405e-04	 &2.0050 & 5.4798e-04 & 1.9981\\
\hline
\end{tabular}
\end{table*}
\\
As shown in Table 2.1, ADI-FDTD is not energy-conserved. It has large errors of energy of $10^{-3} - 10^{-5}$ when compared to the other schemes, EC(1,4), EC(2,4) and EO-EC(1,4). Table 2.1 uses the parameters $\Delta$t = $\Delta$x = $\Delta$y = 0.01 and $k_x = k_y$. As the wavenumber \textbf{k} increases, the errors of energy of EC(1,4), EC(2,4) and EO-4th-EC reaches $10^{-14}$, which is mechanical error. Hence, these three schemes are energy-conserved. \\

Tables 2.2 and 2.3 compare relative accuracy of numerical solutions in time step for the four schemes using parameters, T = 1, $\Delta t =\Delta x =\Delta y=1/N$, $\mu=\epsilon$ = 1, $k_x = k_y$ = 1. It is easily seen that EC(1,4) scheme is only first order in time step, where EC(1,4), ADI-FDTD and EO-4th-EC are all second order in time step. Comparing with the EC(2,4) scheme, our EO-4th-EC scheme takes much less time to run as it only requires two stages at each time step while EC(2,4) uses three stages at each time step.\\

Tables 2.4 and 2.5 compare the relative accuracy in spatial steps for four schemes with parameters, T = 1, $\Delta x =\Delta y=1/N$, $\mu=\epsilon$ = 1, $k_x = k_y$ = 1, $\Delta t=1/N^4$ for EC(1,4) and $\Delta t=1/N^2$ for the other three schemes. The convergent ratio clearly states that ADI-FDTD scheme is second order in space where our scheme EO-4th-EC, and EC(1,4), EC(2,4) schemes are fourth-order accurate in space. Tables 2.2, 2.3, 2.4, and 2.5 prove that the proposed EO-4th-EC scheme is second order in time and fourth order in space.\\

Tables 2.6 and 2.7 list the numerical results of divergence free, Div1 and Div2, for the four schemes mentioned above with $\Delta t=\Delta x=\Delta y$ at time t=1. It can be clearly seen that our EO-4th-EC and EC(2,4), ADI-FDTD schemes are all second order in time while EC(1,4) is first order in time.

\singlespacing
\chapter{Numerical Experiments of Electromagnetic Wave Propagation}
\doublespacing
This chapter shows the numerical experiments of the wave propagation of electromagnetic field in different media by the proposed EO-4th-EC scheme. Wave propagation is analyzed for different cases with Perfectly Electric Conducting (PEC) boundary and with Perfectly Matched Layer (PML) for vacuum domain and for different media in the presence of a source or an initial electric or magnetic field.

\section{Maxwell's Equations with sources in 2D: Model and Scheme}
The EO-4th-EC scheme is used to compute the electromagnetic fields in each domain. For different situations, TE and TM modes can be used. \\
\newpage
\textbf{TE Mode:}\\
\begin{equation}
\frac{\partial{E_x}}{\partial{t}}= \frac{1}{\epsilon} \frac{\partial{H_z}}{\partial{y}} 
\end{equation}
\begin{equation}
\frac{\partial{E_y}}{\partial{t}}= -\frac{1}{\epsilon} \frac{\partial{H_z}}{\partial{x}}
\end{equation}
\begin{equation}
\frac{\partial{H_z}}{\partial{t}}= \frac{1}{\mu}(\frac{\partial{E_x}}{\partial{y}}-\frac{\partial{E_y}}{\partial{x}}) +\frac{1}{\mu} J_{z}.
\end{equation}
\textbf{TM Mode:}\\
\begin{equation}
\frac{\partial{H_x}}{\partial{t}}= \frac{1}{\mu} \frac{\partial{E_z}}{\partial{y}}
\end{equation}
\begin{equation}
\frac{\partial{H_y}}{\partial{t}}= -\frac{1}{\mu} \frac{\partial{E_z}}{\partial{y}}
\end{equation}
\begin{equation}
\frac{\partial{E_z}}{\partial{t}}= \frac{1}{\epsilon}(\frac{\partial{H_x}}{\partial{y}}-\frac{\partial{H_y}}{\partial{x}} )+\frac{1}{\mu} J_{z}.
\end{equation}


In each mode, $J_z$, source in z direction is added on to the $E_z$ or $H_z$ equation. Similarly, if there is a source present in x or y-direction, $J_x$ or $J_y$ term can be added to their corresponding equations.\\$ $\\
Here, $\epsilon$ and $\mu$ are the permittivity and permeability, $E=(E_x(x,y,t)), E_y(x,y,t))$ is the electric field and $H_z=H_z(x,y,t)$ is the magnetic field for $(x,y) \in \Omega$ $ =(0,a)\times (0,b1)$ and $t \in (0,T]$.
Perfectly conducting boundary condition is satisfied on the boundary:
$(E,0) \times (\vec{n},0)$ = 0 on (0,T] $\times \partial{\Omega},$
where $\partial{\Omega}$ is the boundary of $\Omega$ and $\vec{n}$ is the outward normal vector on $\partial{\Omega}$.\\$ $\\
We can split the Maxwell's equations in each time interval $[t^n, t^{n+1}]$ as following:\\
\textbf{TE Mode:}

\begin{equation}
 \frac{\partial{E_x}}{\partial{t}}=\frac{1}{\epsilon}\frac{\partial{H_z}}{\partial{y}},  \nonumber \\
 \end{equation}
\begin{equation} \frac{1}{2}\frac{\partial{H_z}}{\partial{t}}=\frac{1}{\mu}\frac{\partial{E_x}}{\partial{y}} +\frac{1}{2\mu} J_{z}. \\ 
\end{equation} 

\begin{equation}
 \frac{\partial{E_y}}{\partial{t}}=-\frac{1}{\epsilon}\frac{\partial{H_z}}{\partial{x}},  \nonumber \\
\end{equation}
\begin{equation} \frac{1}{2}\frac{\partial{H_z}}{\partial{t}}=-\frac{1}{\mu}\frac{\partial{E_y}}{\partial{x}}+\frac{1}{2\mu} J_{z} . \\ 
\end{equation} 

\textbf{TM Mode:\\}
\begin{equation}
\frac{\partial{H_x}}{\partial{t}}=\frac{1}{\mu}\frac{\partial{E_z}}{\partial{y}},\nonumber
\end{equation}
\begin{equation}
\frac{1}{2}\frac{\partial{E_z}}{\partial{t}}=\frac{1}{\epsilon}\frac{\partial{H_x}}{\partial{y}}+\frac{1}{2\epsilon} J_{z}. 
\end{equation}
\begin{equation}
\frac{\partial{H_y}}{\partial{t}}=-\frac{1}{\mu}\frac{\partial{E_z}}{\partial{x}},\nonumber
\end{equation}
\begin{equation}
\frac{1}{2}\frac{\partial{E_z}}{\partial{t}}=-\frac{1}{\epsilon}\frac{\partial{H_y}}{\partial{x}}+\frac{1}{2\epsilon} J_{z}. 
\end{equation}
The TE mode with source is discretized using the EO-4th-EC Scheme as following.\\

\textbf{Stage 1: The Odd time step.}\\

\textbf{At every odd time step, from $t^{(2k)}$ to $ t^{(2k+1)}$, use $E^{(2k)}_x$, $E^{(2k)}_y$ and $H^{(2k)}_z$ to compute $E^{(2k+1)}_x$, $E^{(2k+1)}_y$ and $H^{(2k+1)}_z$ by the fourth-order spatial scheme.}\\
\textit{Stage 1.1: First, compute $E^{(2k+1)}_y$ and intermediate variable $H^*_z$ using $H^{(2k)}_z$ and $E^{(2k)}_y$; for j=0,1,..,J-1}\\
For the interior points, i=2,3,..,I-2
\begin{multline}
\frac{E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}-E^{(2k)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t}
= - \frac{1}{2\epsilon}\varLambda_x\{H^*_{z_{i,j+\frac{1}{2}}} + H^{(2k)}_{z_{i,j+\frac{1}{2}}}\}, \\
\frac{H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t}
= - \frac{1}{2\mu}\varLambda_x\{E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\} +\frac{1}{2\mu}J^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}},
\end{multline}
For boundary nodes, i=1 and i=I-1,
\begin{multline}
\frac{E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}-E^{(2k)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t}
= - \frac{1}{2\epsilon}\widetilde{\varLambda_x}\{H^*_{z_{i,j+\frac{1}{2}}} + H^{(2k)}_{z_{i,j+\frac{1}{2}}}\}, \\
\frac{H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t}
= - \frac{1}{2\mu}\widetilde{\varLambda_x}
\{E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\}+\frac{1}{2\mu}J^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}},
\end{multline}
\textit{Stage 1.2: Then compute $E^{(2k+1)}_x$ and $H^{(2k+1)}_z$ using $H^{(*)}_z$ and $E^{(2k)}_x$, for i = 0,1,..,I-1}\\
For interior nodes, j=2,3,..,J-2,
\begin{multline}
\frac{E^{(2k+1)}_{x_{i+\frac{1}{2},j}}-E^{(2k)}_{x_{i+\frac{1}{2},j}}}{\varDelta t}
=  \frac{1}{2\epsilon}\varLambda_y\{H^{(2k+1)}_{z_{i+\frac{1}{2},j}} + H^{*}_{z_{i+\frac{1}{2},j}}\} ,\\
\frac{H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t}  =\frac{1}{2\mu}\varLambda_y\{E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\}+\frac{1}{2\mu}J^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}},
\end{multline}
For boundary nodes,j=1 and j=J-1,
\begin{multline}
\frac{E^{(2k+1)}_{x_{i+\frac{1}{2},j}}-E^{(2k)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} = \frac{1}{2\epsilon}\widetilde{\varLambda_y}\{H^{(2k+1)}_{z_{i+\frac{1}{2},j}} + H^{*}_{z_{i+\frac{1}{2},j}}\},\\
\frac{H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{*}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu}\widetilde{\varLambda_y}\{E^{(2k+1)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k)}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}\}+\frac{1}{2\mu}J^{(2k)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}},
\end{multline}
\\
\textbf{Stage 2. The Even time step.}\\

\textbf{At every even time step, from $t^{(2k+ 1)}$ to $t^{(2k+2)}$, use $E^{(2k+1)}_x$, $E^{(2k+1)}_y$ and $H^{(2k+1)}_z$ to compute $E^{(2k+2)}_x$, $E^{(2k+2)}_y$ and $H^{(2k+2)}_z$ by firstly fourth order difference y-direction scheme and secondly fourth order difference x-direction scheme.} \\
\textit{Stage 2.1: First compute $E^{(2k+2)}_x$ and intermediate variable $H^{**}_z$ using $H^{(2k+1)}_z$ and $E^{(2k+1)}_x$, i=0,1,..,I-1.} \\
For interior nodes, j = 2,3,..,J-2
\begin{multline}
\frac{E^{(2k+2)}_{x_{i+\frac{1}{2},j}}-E^{(2k+1)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} =  \frac{1}{2\epsilon}\varLambda_y\{H^{**}_{z_{i+\frac{1}{2},j}} + H^{(2k+1)}_{z_{i+\frac{1}{2},j}}\}, \\
\frac{H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu}\varLambda_y\{E^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\}+\frac{1}{2\mu}J^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}},
\end{multline}
For boundary nodes, j=1 and J-1.
\begin{multline}
\frac{E^{(2k+2)}_{x_{i+\frac{1}{2},j}}-E^{(2k+1)}_{x_{i+\frac{1}{2},j}}}{\varDelta t} =  \frac{1}{2\epsilon}\widetilde{\varLambda_y}\{H^{**}_{z_{i+\frac{1}{2},j}} + H^{(2k+1)}_{z_{i+\frac{1}{2},j}}\},\\
\frac{H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} = \frac{1}{2\mu}\widetilde{\varLambda_y}\{E^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\}+\frac{1}{2\mu}J^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}},
\end{multline}
\textit{Stage 2.2: Then compute $E^{(2k+2)}_y$ and $H^{(2k+2)}_z$ using $H^{**}_z$ and $E^{(2k+1)}_y$, j=0,1,..,J-1}\\
For interior nodes, i=2,3,..,I-2.
\begin{multline}
\frac{E^{(2k+2)}_{y_{i,j+\frac{1}{2}}}-E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t} =  -\frac{1}{2\epsilon}\varLambda_x\{H^{(2k+2)}_{z_{i,j+\frac{1}{2}}} + H^{**}_{z_{i,j+\frac{1}{2}}}\}, \\
\frac{H^{(2k+2)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} =-\frac{1}{2\mu}\varLambda_x\{E^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\}+\frac{1}{2\mu}J^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}},
\end{multline}
For boundary nodes, i=1, and i = I-1:
\begin{multline}
\frac{E^{(2k+2)}_{y_{i,j+\frac{1}{2}}}-E^{(2k+1)}_{y_{i,j+\frac{1}{2}}}}{\varDelta t} =  -\frac{1}{2\epsilon}\widetilde{\varLambda_x}\{H^{(2k+2)}_{z_{i,j+\frac{1}{2}}} + H^{**}_{z_{i,j+\frac{1}{2}}}\},\\
\frac{H^{(2k+2)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}-H^{**}_{z_{i+\frac{1}{2},j+\frac{1}{2}}}}{\varDelta t} =-\frac{1}{2\mu}\widetilde{\varLambda_x}\{E^{(2k+2)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}} + E^{(2k+1)}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}\}+\frac{1}{2\mu}J^{(2k+1)}_{z_{i+\frac{1}{2},j+\frac{1}{2}}},
\end{multline}

The boundary conditions are given by
\begin{equation}
E^*_{x_{i+\frac{1}{2},0}}=E^*_{x_{i+\frac{1}{2},J}}=E^{(2k+2)}_{x_{i+\frac{1}{2},0}}=E^{(2k+2)}_{x_{i+\frac{1}{2},J}}=E^{(2k+2)}_{y_{0,j+\frac{1}{2}}}=E^{(n+1)}_{y_{I,j+\frac{1}{2}}} = 0,
\end{equation}

\newpage

In a case of a dielectric material interface, averages of the media electric permittivitties have been proposed for accurate updating of the EM field components of the interface. Typically, the calculation of the effective permittivities is based on a Taylor series analysis of the electromagnetic field quantities at the dielectric interface. Calculation of effective permittivities for second order accurate finite difference schemes at a dielectric interface is effected in a systematic and straight forward manner through the discretization of the integral forms of Maxwell's equations over finite volumes containing the interface. The proposed analysis in \cite{hwang2001} confirms the validity of the second-order accuracy in time  of the FDTD scheme. The electric permittivity at the interface is defined as:
\begin{equation}
\epsilon=\frac{\epsilon_1 + \epsilon_2}{2}
\end{equation}
\\
In the following section, we run some experimental results of a wave propagation to check the validity of EO-4th-EC scheme in a vacuum and multimedia layer domain. 


\section{Experimental results}
\singlespacing
\subsection{Experiment 1: Wave propagating through vacuum space with Perfectly Electric Conductive Boundary} 
\doublespacing
\subsubsection{Wave propagating through vacuum space with the presence of an initial magnetic field, $H_z_0$}
First we consider a vacuum domain with an initial magnetic field as shown in Figure 3.1. The domain has symmetry about the z-axis and it makes the $\frac{\partial}{\partial z}$ components of the 3D Maxwell's equations equals to zero. We consider initial magnetic field in the z-direction, hence we use the TE model with $E_x$, $E_y$ and $H_z$. First we consider the Perfect Electric Conductive (PEC) Boundary Condition. Later, Perfectly Matched Layer (PML) condition is introduced at the boundary to model the vacuum and multilayer domain.\\

In this experiment, we consider an initial magnetic field in the z-direction, $H_z_0$. The initial magnetic field is stated below and is positioned at the center of the perfectly electric conductive domain. The boundary condition is $E_x_0 = 0, E_y_0 =0$ and $H_z_0$. The source is defined as 
$ \vec{\textbf{J}}=(J_x,J_y, J_z)=0, $ \\ $ $\\
The initial condition is:

\begin{equation}
H_{z_0} = \left\{
\begin{array}{lr}
\frac{1}{25}(10+15 \sin(\frac{10}{3}\pi r)) & r< 0.5 \\
0 & \text{otherwise} 
\end{array}
\right.
\end{equation}
 
\begin{figure}
	\begin{center}
		\begin{tikzpicture}[aspect=0.5]
			\node [cylinder, cylinder uses custom fill, cylinder end fill=red!50,
			cylinder body fill=red!25, rotate=90, minimum height=4.5cm] at (1.5,1.5,1.5) {$H_z_0$};
		\draw (0,0,0) -- (4,0,0)--(4,4,0) --(0,4,0)--(0,0,0)--(0,0,4)--(4,0,4)--(4,0,0)--(0,0,0)--(0,4,0)--(0,4,4)--(0,0,4)--(0,0,0)--(0,4,0)--(4,4,0)--(4,4,4)--(0,4,4)--(0,0,4)--(4,0,4)--(4,4,4);
		\draw[thick,->] (5,0,0) -- (6,0,0) node[anchor=north east]{\textbf{x}};
		\draw[thick,->] (0,4,0) -- (0,5,0) node[anchor=north]{\textbf{z}};
		\draw[thick,->] (0,0,5) -- (0,0,6) node[anchor=south]{\textbf{y}};
		\filldraw [fill=green, draw=black, dashed, fill opacity=0.2] (0,2,0) -- (4,2,0)--(4,2,4)--(0,2,4)--cycle;
		\end{tikzpicture}
		\begin{tikzpicture}
		\filldraw[fill=black, draw=black] (2,2) circle (0.2cm);
		\draw[step=0.5cm,gray,very thin] (0,0) grid (4,4);
		\draw[thick,->] (5,0) -- (6,0) node[anchor= east]{\textbf{x}};
		\draw[thick,->] (5,0) -- (5,1) node[anchor= north]{\textbf{y}};
		\end{tikzpicture}
		\caption{(Experiment 1) Initial magnetic field, $H_z_0$, in the center of vacuum 3D (top) and cross-section of x-y plane(bottom).}
	\end{center}
\end{figure}
Figure 3.1 includes the domain of 5.08mm $\times$ 5.08mm $\times$ 5.08mm  vacuum domain. The dielectric constant, i.e., relative electric permittivity and the relative magnetic permeability is 1. Figure 3.2 shows the mesh figures of magnetic field in the z-direction as it travels from T = 0 to T=4. As time approaches 4, the magnetic wave hits the perfectly electric conductive boundary, and it bounces back producing a reflection. The result behaves as it was expected and stated in \cite{Inan2011}. Figure 3.3 shows the contour figures of electric field in x-direction, y-direction and the magnetic field in z-direction as it travels to T=3.14. The $H_z$ figure clearly slows the reflection at each side of the domain.
\begin{figure}
	\includegraphics[width=0.5\textwidth]{1a}
	\includegraphics[width=0.5\textwidth]{1T=1}\\
	\includegraphics[width=0.5\textwidth]{1T=2}
	\includegraphics[width=0.5\textwidth]{1T=4}\\
	\caption{Mesh figures of magnetic field,$H_z$ in the z-direction (Top) Left: T=0, Right: T=1, (Bottom) Left: T=2, Right: T=4.}
\end{figure}

\begin{figure}
	\begin{center}
\includegraphics[width=0.63\textwidth]{PECSameMedia_Ex}
\includegraphics[width=0.63\textwidth]{PECSameMedia_Ey}\\
\includegraphics[width=0.63\textwidth]{PECSameMedia_Hz}\\
	\caption{Contour Figures for Ex, Ey, Hz at T=3.14.}
\end{center}
\end{figure}	
\newpage
\singlespacing
\subsection{Experiment 2: Wave propagating through different media with Perfectly Electric Conductive Boundary}

\doublespacing
\subsubsection{Case a: Wave propagating through different media with a higher electric permittivity in the center}

Figure 3.4 displays the domain of an initial magnetic field as it travels through vacuum to a conductive medium. The electric permittivitty and magnetic permeability is modified at the interface to have second order in time. We take the arithmetic mean at the interface for relative electric permittivity as defined in equation (3.24) as shown in paper \cite{Hirono2000}. The relative magnetic permeability is defined  similarly. This experiment also carries a symmetry in the z-direction, hence TE mode is used with Electric field in the x and y-direction and Magnetic field in the z-direction, $E_x$, $E_y$ and $H_z$. PEC boundary condition is considered. We use the same initial magnetic field as used in the previous experiment as stated in equation (3.21).\\
\begin{figure}
	\begin{center}
		\begin{tikzpicture}		
		\draw[black,thick] (0,0,0) -- (4,0,0)--(4,4,0) --(0,4,0)--(0,0,0)--(0,0,4)--(4,0,4)--(4,0,0)--(0,0,0)--(0,4,0)--(0,4,4)--(0,0,4)--(0,0,0)--(0,4,0)--(4,4,0)--(4,4,4)--(0,4,4)--(0,0,4)--(4,0,4)--(4,4,4);
		\filldraw[fill=yellow, draw=black] (1,4,0)--(1,4,4)--(0,4,4)--(0,4,0)--cycle;
		\filldraw[fill=yellow, draw=black] (0,4,4)--(1,4,4)--(1,0,4)--(0,0,4)--cycle;
		\filldraw[fill=yellow, draw=black] (1,4,0)--(1,4,4)--(1,0,4)--(1,0,0)--cycle;
		
		\node[cylinder, draw=black,thick, fill=pink, rotate=90, shape aspect=1, minimum height=4cm] at (2,2,2);
		
		\filldraw[fill=yellow, draw=black] (3,4,0)--(3,4,4)--(4,4,4)--(4,4,0)--cycle;
		\filldraw[fill=yellow, draw=black]  (3,0,4)--(4,0,4)--(4,4,4)--(3,4,4)--cycle;
		\filldraw[fill=yellow, draw=black] (4,4,0)--(4,4,4)--(4,0,4)--(4,0,0)--cycle;
				
		\filldraw [fill=green, draw=black, dashed, fill opacity=0.2] (0,2,0) -- (4,2,0)--(4,2,4)--(0,2,4)--cycle;
		\draw[thick,->] (5,0,0) -- (6,0,0) node[anchor=north east]{\textbf{x}};
		\draw[thick,->] (0,4,0) -- (0,5,0) node[anchor=north]{\textbf{z}};
		\draw[thick,->] (0,0,5) -- (0,0,6) node[anchor=south]{\textbf{y}};
		\end{tikzpicture}
		\begin{tikzpicture}
		\draw  (0,0)--(4,0)--(4,4)--(0,4)--cycle;
		\filldraw[fill=yellow, draw=black] (0,0)--(1,0)--(1,4)--(0,4)--cycle;
		\filldraw[fill=yellow, draw=black] (3,0)--(4,0)--(4,4)--(3,4)--cycle;
		\filldraw[fill=pink, draw=black] (2,2) circle (0.2cm);
		\draw[step=0.5cm,gray,very thin] (0,0) grid (4,4);
		\draw[thick,->] (5,0) -- (6,0) node[anchor= east]{\textbf{x}};
		\draw[thick,->] (5,0) -- (5,1) node[anchor= north]{\textbf{y}};
		\end{tikzpicture}
		\caption{(Experiment 2) Initial Magnetic field in the z-direction, testing the behaviour inside different media in 3D (top), and 2D (bottom).}
	\end{center}
\end{figure}
\\
At the interface, electric permittivitty and magnetic permeability is defined as follows:

Relative Electric permittivity:\\
\[ 	\epsilon = \begin{cases} 
1 & a_1 < x < i_1 \\
12 & i_1 < x < i_2 \\
1 & i_2 < x <b_1 \\
\end{cases}
\]
\begin{equation}
\end{equation}\\
By paper \cite{Hirono2000}, at the interface, $\epsilon $ is defined as \\
\[ 	\epsilon = \begin{cases} 
\frac{1+12}{2} &  x = i_1 \\
\frac{1+12}{2} & x= i_2 \\\end{cases}
\]
\begin{equation}
\end{equation}
Here, $i_1$ and $i_2$ are interfaces in Figure 3.4. $\mu$ can be defined similarly at the interface. 

\begin{figure}
	\includegraphics[width=0.5\textwidth]{2a}
	\includegraphics[width=0.5\textwidth]{2T=1}\\
	\includegraphics[width=0.5\textwidth]{2T=2}
	\includegraphics[width=0.5\textwidth]{2T=5}\\
	\caption{Magnetic Field, $H_z$, in different Domains, Top: Left: T=0, Right: T=1, Bottom: Left: T=2, Right: T=5.} \\$ $\\
\end{figure}


\begin{figure}
	\begin{center}
\includegraphics[width=0.7\textwidth]{PEC_Ex_T=4}
\includegraphics[width=0.7\textwidth]{PEC_Ey_T=4}
\includegraphics[width=0.7\textwidth]{PEC_Hz_T=4}
	\caption{Contour Figures for Ex, Ey, Hz at T=4.}
	\end{center}
\end{figure}
Figure 3.5 glances at the mesh figures of magnetic field in the z-direction for T=0 to T=5. The wave travels from a metallic sheet located at the center of the domain to vacuum domain which is located on the left and right side of the domain. Since the relative electric permittivity is higher in the center, it can be clearly noted that the magnetic field is stronger in the center and it decreases as it moves outward to air domain, where the electric permittivity is lower. In the long run, the magnetic field is stronger in the metallic domain. Figure 3.6 provides the contour figures for $E_x$, $E_y$ and $H_z$ at T=4. The contour figures clearly show that the wave changes shape as the wave propagates from metal to air.



\subsubsection{Case b: Wave propagating through different media with a lower electric permittivity in the center}
Further, we introduce a gold sheet on the left and right vacuum at the center of domain. Same domain, symmetry and initial magnetic fields as experiment 2a are considered. We continue to use the TE mode with electric field in the x and y-direction with magnetic field in the z direction, $E_x$, $E_y$ and $H_z$. \\


Here, the relative electric permittivity is $\epsilon_r$= $\epsilon / \epsilon_0$, where $\epsilon_0 = 8.854\times10^{-7}$\\
\[ 	\epsilon_r = \begin{cases} 
6.9 & a1 < x < i1 \qquad \text{(Gold)}\\
1 & i1 < x < i2 \qquad \text{(Vacuum)}\\
6.9 & i2 < x <b1 \qquad \text{(Gold)} \\
\end{cases}
\]
\begin{equation}
\end{equation}
At the interface, $\epsilon_r$ is defined as \\
\[ 	\epsilon_r = \begin{cases} 
\frac{1+6.9}{2} &  x = i1 \\
\frac{1+6.9}{2} & x= i2 \\
\end{cases}
\]
\begin{equation}
\end{equation}

\text{$\mu$ can be defined similarly at the interface.}
\\
\begin{figure}
	\includegraphics[width=0.5\textwidth]{3T=1}
	\includegraphics[width=0.5\textwidth]{3T=3}\\
	\includegraphics[width=0.5\textwidth]{3T=5}\\
	\caption{Mesh Figures of magnetic Field, $H_z$ Passing from vacuum to gold, Top: T=1, Center: T=3, Bottom: T=5.} \\$ $\\
\end{figure}

Figure 3.7 shows the mesh of magnetic field in the z-direction, $H_z$, as the initial magnetic field travels from vacuum to gold plates. With increase in time, the magnetic field is stronger in the gold plates due to higher relative electric permittivity.


\singlespacing
\subsection{Experiment 3: Wave propagating through vacuum domain with a Perfectly Matched Layer}
\doublespacing
In this experiment, we consider the same domain listed in Figure 3.1, but it is now surrounded by a Perfectly Matched Layer. There is a symmetry in the z direction. Hence, TE mode with Electric field in the x, y-direction, $E_x, E_y$, and Magnetic field in the z-direction$H_z$ is used. 

\subsubsection{Case a: Initial Magnetic Field, $H_z0$}
First, an initial magnetic field in the z direction, $H_z_0$, is introduced as shown in Figure 3.8 (Left) and the magnetic field in the z-direction, $H_z$, is reviewed over time.\\
\begin{figure}
		\begin{tikzpicture}
		\filldraw[fill=pink, draw=pink] (-0.5,-0.5)--(-0.5,4.5)--(0,4.5)--(0,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (4,-0.5)--(4,4.5)--(4.5,4.5)--(4.5,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (0,-0.5)--(0,0)--(4,0)--(4,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (0,4)--(0,4.5)--(4,4.5)--(4,4)--cycle;
		\filldraw[fill=black, draw=black] (2,2) circle (0.2cm);
		\draw[step=0.5cm,gray,very thin] (0,0) grid (4,4);
		\draw[thick,->] (5,0) -- (6,0) node[anchor= east]{\textbf{x}};
		\draw[thick,->] (5,0) -- (5,1) node[anchor= north]{\textbf{y}};
		\filldraw[fill=pink, draw=pink] (-0.5+7,-0.5)--(-0.5+7,4.5)--(0+7,4.5)--(0+7,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (4+7,-0.5)--(4+7,4.5)--(4.5+7,4.5)--(4.5+7,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (0+7,-0.5)--(0+7,0)--(4+7,0)--(4+7,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (0+7,4)--(0+7,4.5)--(4+7,4.5)--(4+7,4)--cycle;
		\filldraw[fill=green, draw=black] (2+7,2) circle (0.2cm);
		\draw[step=0.5cm,gray,very thin] (0+7,0) grid (4+7,4);
		\draw[thick,->] (5+7,0) -- (6+7,0) node[anchor= east]{\textbf{x}};
		\draw[thick,->] (5+7,0) -- (5+7,1) node[anchor= north]{\textbf{y}};
		\end{tikzpicture}
\caption{(Experiment 3)(Left) Case a - 2D Domain for a vacuum space with Perfectly Matched Layer with $H_z_0$, (Right) Case b - 2D Domain for a vacuum space with Perfectly Matched Layer with an initial source in the center, $J_z$.}
\end{figure}

\begin{figure}
        \includegraphics[width=0.5\textwidth]{PMLT=1}
		\includegraphics[width=0.5\textwidth]{PMLT=3}\\
		\includegraphics[width=0.5\textwidth]{PMLT=4}
		\caption{{Hz with Initial Magnetic Field}, Top: Left T=1, Right T=3, Bottom: T=4.}
\end{figure}	
\begin{figure}
	\includegraphics[width=0.4\textwidth]{PMLST=25}
	\includegraphics[width=0.4\textwidth]{PMLST=3}\\
	\includegraphics[width=0.4\textwidth]{PMLST=4}\\
	\caption{Hz with initial Source in z-direction in air domain $J_x =0, J_y =0, J_z=A \sin(2\pi f_0t)$,  Top: T=2.5, Center: T=3, Bottom: T=4.}
\end{figure}
\\ 

Figure 3.9 shows the magnetic field in z direction in a case where initial magnetic field, $H_z_0$, is present. As seen in this figure, the magnetic wave starts to travel outwards with increase in time. As the wave hits the perfectly matched layer, the wave smoothly makes an exit by being absorbed into the layer and creating no reflection.


\subsubsection{Case b: Initial Electric Source in z direction, $J_z$}	
Next, we excite our domain, Figure 3.8 (right), with a magnetic current source in the z-direction which is defined as a sinusoidal function. The initial magnetic current source is defined as:
$$\textbf{J}=(\vec{J_x},\vec{J_y} \text{ and } \vec{J_z}) =(0,0,J_z),$$ where
$$ J_z(x_0,y_0) = A \sin(2 \pi f_0 t),$$
here, t is time, A is the amplitude, $f_0$ is the frequency of the wave. The magnetic current source is located at the center of the domain. Figure 3.10 shows the magnetic field in z-direction with a magnetic current source, located at the center. As time increases, the magnetic waves keeps on moving outwards and is absorbed into the perfectly matched layer.

\subsection{Experiment 4: Wave Propagating in different media with a Perfectly Matched Layer}
This section focuses on the behaviour of a traveling EM wave inside different media whose domain is surrounded by a Perfectly Matched Layer. The EO-4th-EC scheme is simulated in MATLAB. Symmetry in the z axis is considered and hence, TE mode with $E_x$, $E_y$ and $H_z$ is used. In this case, we only consider the presence of an initial magnetic field in the z direction. Figure 3.11 shows the cross-section of the x-y plane with the presence of a Perfectly matched layer. Initial magnetic field is defined as equation (3.21).
\subsubsection{ $H_z$ in z-direction Vacuum(center) to Gold (left and right)}
\begin{figure}
	\begin{center}
		\begin{tikzpicture}
		\draw  (0,0)--(3,0)--(3,3)--(0,3)--cycle;
		\filldraw[fill=yellow, draw=black] (0,0)--(1,0)--(1,3)--(0,3)--cycle;
		\filldraw[fill=yellow, draw=black] (2,0)--(3,0)--(3,3)--(2,3)--cycle;
		\filldraw[fill=pink, draw=black] (1.5,1.5) circle (0.2cm);
		\draw[step=0.5cm,gray,very thin] (0,0) grid (3,3);
		\filldraw[fill=pink, draw=pink] (-0.5,-0.5)--(-0.5,3.5)--(0,3.5)--(0,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (3,-0.5)--(3,3.5)--(3.5,3.5)--(3.5,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (0,-0.5)--(0,0)--(3,0)--(3,-0.5)--cycle;
		\filldraw[fill=pink, draw=pink] (0,3)--(0,3.5)--(3,3.5)--(3,3)--cycle;
		\end{tikzpicture}
		\caption{(Experiment 4) 2D Domain of different media with initial magnetic field with a PML layer.}
	\end{center}
\end{figure}
\begin{figure}
	%	\begin{center}
	\includegraphics[width=0.4\textwidth]{PMLDT=1}
	\includegraphics[width=0.4\textwidth]{PMLDT=3}\\
	\includegraphics[width=0.4\textwidth]{PMLDT=4}
	\includegraphics[width=0.4\textwidth]{PMLDT=5}\\
	\caption{{Hz with initial magnetic field, $J_z$ in z-direction Vacuum(center) to Gold(left and right)}Top Left: T=1, Right: T=3, Bottom: Left: T=4, Right: T=5.}
	%	\end{center}
\end{figure}

\begin{figure}
	\includegraphics[width=0.4\textwidth]{PMLReal_Ex}
	\includegraphics[width=0.4\textwidth]{PMLReal_Ey}\\
	\includegraphics[width=0.4\textwidth]{PMLReal_Hz}
	\caption{Contour for Ex, Ey and Hz at T=5.}
\end{figure}
Figure 3.12 shows the magnetic field in z direction as the initial magnetic field travels through the domain. Bigger fraction of the magnetic field moves to the gold plate rather than traveling in the vacuum domain due to its higher electric permittivity. As it hits the boundary, the wave is absorbed into the domain making a smooth exit. Figure 3.13 shows the contour figures of $E_x$, $E_y$ and $H_z$ at T=5. It can be easily seen that the wave changes as it moves to a different domain. This figure also shows the smooth exit each wave makes as it reaches the boundary. 
\subsubsection{PEC versus PML condition}
\begin{figure}
	\begin{center}
	\includegraphics[width=0.7\textwidth]{PEC_Hz_T=5}
	\includegraphics[width=0.7\textwidth]{PML_Hz_T=5}\\
	\caption{PEC vs PML boundary condition for Hz at T=5.}
	\end{center}
\end{figure}
Figure 3.14 compares the magnetic field in the z-direction, $H_z$, at T=5 for a PEC (left) and PML (right) boundary condition. For the PEC boundary domain, the magnetic field bounces the boundary walls and reflects back creating a EM field inside the domain. When comparing it to the PML domain, when the magnetic field reaches the boundary, the EM waves are absorbed in to the PML domain and there is little/no reflection present.


\newpage
\singlespacing
\chapter{Modeling of Coplanar Waveguides}
\doublespacing
Maxwell's equations are commonly used to compute the electric and magnetic fields in different media. They are used to model transmission lines such as microstrip, rectangular waveguide, striplines, coplanar waveguides and etc. They are also used to transmit waves for communication purposes. In this chapter, we will focus on coplanar waveguides and their varying shapes.
Section 4.2 focuses on the method of lumped ports to excite a coplanar waveguide. The domain of a coplanar waveguide is considered to 4 different simulations to see the behaviours of the propagating waves in different media. Further, a cylindrical shaped source is introduced inside a grounded coplanar waveguide. We analyze the benefits and the wave propagation inside the elevated coplanar waveguide. Furthermore, we consider and analyzed a complex transition between a coplanar waveguide and a rectangular waveguide where sources are introduced at different locations.

\section{Coplanar Waveguide Problems}
 Waveguide is a structure which guides electromagnetic waves. Most common waveguides contain metal pipe/line used to carry high frequencies, radio waves particularly microwaves. Metallic lines are used at microwave frequencies utilized to connect receivers and transmitters with antennas. The walls of a waveguide are made of conducting material not to conduct energy but for reflection of waves. Transmission lines are commonly used with sampling circuits for millimeter-wave instruments. CPWs are used in integrated circuits. Integrated circuits are used in all electronic equipment such as computers, mobile phones and other digital home appliances.  CPWs are considered for use for millimeter-wave integrated circuits(MMIC) instead of microstrip. MMIC is an integrated circuit that operates at microwave frequencies of 300MHz to 300GHz. Inputs and Outputs are matched to Characterisitic Impedance (Z)=50 Ohms in a microwave frequency range. Coplanar waveguides are fabricated using semiconductor, Gallium Arsenide(GaAS). Other materials used are Si (Silicon), InP (Indium Phosphide), SiGe (Silicon Germanium), and GaN (Gallium Nitride). Depending on the output required, different fabrication is used. GaAs is the most common dielectric substrate for most transmission lines.\\

\begin{figure}	
		\includegraphics[width=0.5\textwidth]{integrated}
		\includegraphics[width=0.5\textwidth]{2aa}\\
		\caption{Integrated Circuit \& a CPW.}
\end{figure}

Coplanar waveguide was initially introduced by Cheng P Wen in 1969. He first decided to name it ``Planar Stripline", however his coworker suggested him to name it coplanar waveguide, CPW, after his initials. Coplanar waveguide, CPW, has two parallel plates made of conducting material, such as copper or gold, that run with some dielectric materials in between. In an ideal case, the thickness of a dielectric material is infinite. However, in practice, the dielectric substrate is thick enough that the electromagnetic waves die out before they leave the dielectric material. A dielectric medium can sustain electric field but does not conduct electricity. The advantage of a coplanar waveguide is that active devices can be mounted on top of a CPW. It can also provide a very high frequency response.  Coplanar Waveguide can be excited in few ways. The electromagnetic waves transmit inside a waveguide either coming from a source outside or from a source implanted inside the waveguide. A lumped port can be attached to one side of the waveguide to its conducting plates to excite the waveguide. In some cases, a coaxial cable is also attached to the coplanar waveguide depending on the type of situation a coplanar waveguide is used for. Electromagnetic field can be transmitted from antennas to waveguides and can also be transmitted through waveguides.  Paper \cite{Jackson1986}  suggested that CPW is an alternative to microstrip since it works better at millimeter-wave frequencies. CPW was not used due to mistaken assumption that it has a higher conduction loss than microstrip and via holes are not necessary and conductors do not have to be made excessively thin. However, conduction loss depends on conductor's surface roughness and the behaviour of current crowded near the edge of the conductor. \cite{Jackson1986} paper compared the conductor loss for CPW and microstrip in terms of full wave analysis. The paper also optimized the size of CPW in order to compete with microstrip in terms of loss over a large impedance range of a semiconductor.\\

\begin{figure}
	\includegraphics[width=0.5\textwidth]{CPW}			
	\includegraphics[width=0.5\textwidth]{microstrip}\\
	\caption{ CPW (Left), Microstrip (Right).}
\end{figure}

Full Wave Analysis has been applied to find the loss/dispersion analysis of CPW. The analysis was further used to find the impedance and work out the discontinuities at variables or wherever the system is transiting from one device to another. The analysis has been used to find propagation constant, stability condition and used to compare frequency with other parameters such as dielectric constant and phase velocity. Coplanar Waveguides, just like other transmission line, has radiation, dielectric and conduction losses. Lower the overall loss, better a transmission line is considered. Each characteristic Impedance correspond to different dimensions of CPW along with a certain frequency. In paper \cite{Jackson1986}, full wave analysis was used to compare the conductor loss, dispersion and radiation into parasitic modes for a CPW transmission line and a microstrip. The paper showed that the dimensions of CPW are to be chosen to give better results in terms of loss and high frequency response. Although microstrip is well understood at microwave frequencies and is flexible so that large number of circuit elements can be made with it, but they are not the best choice for integrated circuits operating at millimeter wave frequencies. Their disadvantage is that via holes are required to ground active devices. At millimeter wave frequencies, these vias can create significant inductance and degrade circuit performance. CPW is well suited for use with field-effect transistors, especially at millimeter wave frequencies where RF grounding must be close to the device. Ground connections can conveniently be made at the substrate edge. The paper showed that conductor loss depends on conductor surface roughness which varies from substrate to substrate. It also depends upon the behaviour of current crowded near the edge of etched conductor. In paper \cite{Guo2004}, Leakage loss of a mode was represented by the imaginary part of its complex propagation constant. To model the leakage, open boundaries condition is used and it produces no reflection at the boundary. PML Boundary Condition is mostly efficient; However in order to compute the leakage loss from the material loss/gain, PML layer is removed and PEC boundary condition is used. Real dielectric constant is replaced with a complex one whose imaginary part represent the loss/gain of the media. The calculated propagation constant reflects the material loss.\\


In this chapter, we will implement a current source inside the coplanar waveguide and see how the electromagnetic field changes as it moves in different media. Later, we can discover how the magnetic field is transmitted from a coplanar waveguide to a rectangular waveguide. A current source is added in the coplanar waveguide at the point $(x_0,y_0)$ which moves in a certain direction depending on the domain. A sinusoidal wave is used to excite the coplanar waveguide plane at a frequency  in the LMDS frequency range. The source term can be defined as $J_k = \sin(2\pi f_{0}n\Delta t)$, where k=x,y or z. Here, n is the number of steps used in the mesh  and $\Delta t = \frac{T}{n}$, where T is the time for the simulation is run for and $f_0$ is the frequency of the sinusoidal wave traveling inside the waveguide.\\

	
	\begin{figure}
		\includegraphics[width=0.5\textwidth]{View}
		\includegraphics[width=0.5\textwidth]{geometry}\\
		\caption{(Left) Schematic View of the proposed transition of CPW to rectangular waveguide; (Right) Geometry of the CPW-Rectangular waveguide transition, w=0.762mm, g=0.254mm, l=1.092mm, L=2.413mm, s=0.249mm, d=0.191mm, a=5.08mm.}		
	\end{figure}

 Later we can discover the change in the electromagnetic field as the wave propagates from a coplanar waveguide to a rectangular waveguide. A transition working in the LMDS frequency range (26 GHz to 29 GHz) is considered. The source can be in the form of current, voltage implanted or some initial electromagnetic wave from outside the domain.  The EO-4th-EC scheme can be used to compute the waves traveling inside the waveguide as the source creates the electromagnetic field. In the past, microstrips were commonly used in printed circuit boards. CPWs are alternatives to microstrip lines and are super conductive but not as commonly used due to lack of research done on CPWs. Paper \cite{Deslandes2001} looked at the integration of active and passive components through rectangular waveguides which usually requires transitions from planar to non-planar circuits. Figure 4.3 shows the schematic view of the proposed transition in paper \cite{Deslandes2001}. Various approaches had been made for such transition but they all led to complex mounting structures.  If a microstrip circuit is considered in the transition, it is often needed to be cut into a specific shape which is hard to achieve in millimeter-wave range. If a Rectangular waveguide is considered, its components are long and expensive to manufacture and they make the planar/nonplanar integration bulky and costly. Hence, an idea of composing a CPW and rectangular waveguide was proposed. In the paper, back to back transition structure is required from CPW to rectangular waveguide. The waveguide is synthesized with a linear array of metalized via holes on the same substrate used for planar circuit. Via holes are added to remove potential parallel plate mode that can propagate between CPW and the back plane before the rectangular Waveguide. The waveguide can be integrated with metalized wall which allows a very compact integration. An efficient transition between both structures was required hence a stub is added to connect the two waveguides efficiently. This type of system is used in high speed electronic circuits. The paper further compared the power loss vs frequency as the wave travels from CPW to rectangular waveguide.\\


In general, coplanar waveguides are considered due to its property of having low loss comparing to microstrip. Losses are combined of dielectric loss, conduction loss and radiation loss. Figure 4.2 looks at the schematic and geometric view of the proposed transition in paper \cite{Deslandes2001} of coplanar waveguide and rectangular waveguide integrated on the same surface. This example will be solved in Section 4.5. Figure 4.3 looks at the geometry of a CPW and compares to a microstrip. Coplanar waveguide consists of parallel ground plates that are made of conducting material where a microstrip only has one conducting ground plane. Both are placed on top of a dielectric substrate. Due to geometry of CPW, it has lower loss than a microstrip. Nonlinear Transmission lines are commonly used together with sampling circuits for millimeter-wave instruments. In paper \cite{Hofschen1996}, elevated coplanar waveguide was proposed to increase with bandwidth of sampling circuits. Increasing the elevation height of 3$\mu$m decreased the overall loss by almost half for the whole frequency range. The paper also mentioned that phase velocity is approximately doubled. The geometry of an elevated coplanar waveguide is found in Figure 4.6.\\

	\begin{figure}
\includegraphics[width=1.2\textwidth]{MicrostripLoss}\\
\caption{Conductor (Left) and Dielectric (Right) Loss in a Microstrip.}
\includegraphics[width=1.2\textwidth]{CPWLoss}\\
	
\caption{Conductor (Left) and Dielectric (Right) Loss in a CPW.}

	\end{figure}
	\begin{figure}
		\begin{center}
			\includegraphics[width=0.7\textwidth]{example2}\\
			\caption{ Elevated CPW.}
		\end{center}
	\end{figure}
	
	
	
	\begin{figure}
		\includegraphics[width=1.2\textwidth]{example1loss}\\
		\caption{Dielectric (Left) and Conductor (Right) Loss for CPW vs Frequency of a nonelevated CPW.}
		\includegraphics[width=1.2\textwidth]{example2loss}\\
		\caption{ Dielectric (Left) and Conductor (Right) Loss for elevated CPW vs Frequency when the center conductor is elevated by 3mm.}
	\end{figure}

Dielectric Loss is an electric energy converted into heat in a dielectric system subjected to a varying electric field. Conductor Loss is loss occurring in a conductor due to the flow of current.In \cite{Jackson1986}, Full wave analysis was used to compare the conductor loss, dispersion and radiation into parasitic modes for a coplanar waveguide. Microstrip is well understood at microwave frequencies and is flexible so that large number of circuits are made with it. However, they are not the best choice for high frequencies. CPW is an alternate to microstrip but it is not widely used due to the mistaken assumption that it has a higher conduction loss. However, due to research, it was found that the conduction loss solely depends on the roughness of the conductor surface. If CPW and a microstrip are both made of the same conducting material with the same roughness, CPW has a overall lower loss, hence, a better transmission line. Most power losses will be noted as an increase in heat in the specified material. Larger diameter conductors will provide less conductor loss. To ensure the concept of loss behind these papers, losses were computed in microstrip and coplanar waveguides. As shown in the Figures 4.4 and 4.5, CPW has a lower dielectric loss comparing to Microstrip. However CPW and Microstrip has the same conductivity loss if they are made of the same material. Overall, if the two losses are added CPW provides us with less loss and hence is better to use due to its geometry. Looking at Figures 4.7 and 4.8, you notice that the dielectric Loss for coplanar waveguide is -0.0137 but the dielectric Loss for Elevated CPW = $-1.0152 \times 10^{-5}$. The elevated coplanar waveguide gives us a very low dielectric loss comparing to a regular coplanar waveguide.\\




\section{Current due to two lumped ports}
\textbf{Different ways to Excite a Coplanar Waveguide}\\
There are different ways to excite a CPW. Adding a rectangular faces to the model, either normal or parallel to the plane of the CPW, along with a lumped port is one of the way to excite a CPW. A lumped port excitation is applied between the rectangular face and the center conductor. The arrows apply to show the direction of the flow of applied current that behaves sinusoidally dependent of time as seen in the top two figures of Figure 4.9. Another simpler method is to include two lumped ports which either both travel towards or away from the center conductor. In this scenario, the structure is simpler, however two lumped ports are required to be manually set. Figure 4.9 shows each of these methods of excitation. A lumped port is a transmission connection that applies uniform electric field between two metallic boundaries. The excitation at the port can be expressed as a voltage or a current. In this section, we will use electric current as a medium of excitation.\\
\begin{figure}
	\includegraphics[width=1\textwidth]{1aa}\\
	\includegraphics[width=0.5\textwidth]{2aa}
	\caption{Different Methods to excite a CPW: Top: Rectangular face(Normal or Parallel), Center: Two Lumped ports.}
\end{figure}



Now, we focus on running experiments to see the behaviour of an electric current. We start with a simple domain where, only one metallic plate with two dielectric slabs on the sides are connected. Further, the other aspects of the coplanar waveguide are added and analyzed. Each system has a symmetry in the y-direction, hence, we set $\frac{\partial}{\partial y}$ terms to 0. Simplifying, we use the TE model with electric field in the x and z direction, $E_x$, $E_z$ and a magnetic field in the y-direction, $H_y$.\\


\begin{figure} {\textit{\textbf{Simulation 1 Domain}}}\\
\begin{tikzpicture}
\filldraw[fill=pink, draw=black] (1,0.3,0)--(1,0.3,4)--(0,0.3,4)--(0,0.3,0)--cycle;
\filldraw[fill=pink, draw=pink] (0,0.3,4)--(1,0.3,4)--(1,0,4)--(0,0,4)--cycle;
\filldraw[fill=pink, draw=pink] (1,0.3,0)--(1,0.3,4)--(1,0,4)--(1,0,0)--cycle;

\filldraw[fill=cyan, draw=cyan] (1,0,4)--(3,0,4)--(3,0.3,4)--(1,0.3,4)--cycle;
\filldraw[fill=cyan, draw=cyan] (1,0.3,4)--(3,0.3,4)--(3,0.3,0)--(1,0.3,0)--cycle;

\filldraw[fill=pink, draw=pink] (3,0.3,0)--(3,0.3,4)--(4,0.3,4)--(4,0.3,0)--cycle;
\filldraw[fill=pink, draw=pink]  (3,0,4)--(4,0,4)--(4,0.3,4)--(3,0.3,4)--cycle;
\filldraw[fill=pink, draw=pink] (4,0.3,0)--(4,0.3,4)--(4,0,4)--(4,0,0)--cycle;

\draw[thick,->] (4,0,0) -- (5,0,0) node[anchor=north east]{\textbf{x}};
\draw[thick,->] (0,0.5,0) -- (0,1,0) node[anchor=north]{\textbf{z}};
\draw[thick,->] (0,0,5) -- (0,0,6) node[anchor=south]{\textbf{y}};

\draw [red,thick,->] (0,0.3,0) -- (1,0.3,0);
\draw [red,thick,->] (0,0.3,0.5) -- (1,0.3,0.5);
\draw [red,thick,->] (0,0.3,1) -- (1,0.3,1);
\draw [red,thick,->] (0,0.3,1.5) -- (1,0.3,1.5);
\draw [red,thick,->] (0,0.3,2) -- (1,0.3,2);
\draw [red,thick,->] (0,0.3,2.5) -- (1,0.3,2.5);
\draw [red,thick,->] (0,0.3,3) -- (1,0.3,3);
\draw [red,thick,->] (0,0.3,3.5) -- (1,0.3,3.5);
\draw [red,thick,->] (0,0.3,4) -- (1,0.3,4);

\draw [red,thick,<-] (3,0.3,0) -- (4,0.3,0);
\draw [red,thick,<-] (3,0.3,0.5) -- (4,0.3,0.5);
\draw [red,thick,<-] (3,0.3,1) -- (4,0.3,1);
\draw [red,thick,<-] (3,0.3,1.5) -- (4,0.3,1.5);
\draw [red,thick,<-] (3,0.3,2) -- (4,0.3,2);
\draw [red,thick,<-] (3,0.3,2.5) -- (4,0.3,2.5);
\draw [red,thick,<-] (3,0.3,3) -- (4,0.3,3);
\draw [red,thick,<-] (3,0.3,3.5) -- (4,0.3,3.5);
\draw [red,thick,<-] (3,0.3,4) -- (4,0.3,4);

	\filldraw [fill={rgb:black,1;white,2}, draw=black, dashed, fill opacity=0.2] (0,0.3,2) -- (4,0.3,2)--(4,0,2)--(0,0,2)--cycle;


\draw  (0+7,0-2)--(4+7,0-2)--(4+7,1.5-2)--(0+7,1.5-2)--cycle;
\filldraw[fill=pink, draw=pink] (0+7,0-2)--(1+7,0-2)--(1+7,1-2)--(0+7,1-2)--cycle;
\filldraw[fill=cyan, draw=cyan] (1+7,0-2)--(3+7,0-2)--(3+7,1-2)--(1+7,1-2)--cycle;
\filldraw[fill=pink, draw=pink] (3+7,0-2)--(4+7,0-2)--(4+7,1-2)--(3+7,1-2)--cycle;
\draw[step=0.5cm,gray,very thin] (0+7,0-2) grid (4+7,1.5-2);
\draw [red, thick,->] (0+7,1-2) -- (1+7,1-2);
\draw [red, thick,<-] (3+7,1-2) -- (4+7,1-2);
\draw[thick,->] (4+7,0-2) -- (4.5+7,0-2) node[anchor=east]{\textbf{x}};
\draw[thick,->] (0+7,1.5-2) -- (0+7,2.5-2) node[anchor=north]{\textbf{z}};
\end{tikzpicture}

\caption{Left: Metallic Slab with two Dielectric slabs attached along with current running in the x direction, Right: 2D Domain of the above. }
\end{figure}

Figure 4.10 is the simulation domain of one metallic plate surrounded by two parallel dielectric slabs with current carrying lumped ports. The lumped ports can carry electric current either towards or away from the metallic plate. In this example, we first run electric current of a line source that is only time dependent. 
$$J_x = \sin(2\pi f_0 t)$$
Then, we further run a line source which is both time and space dependent.
$$J_x = \sin(2\pi f_0 x) \times \sin(2\pi f_0 t)$$

\begin{figure}
\includegraphics[width=0.55\textwidth]{1LineEz_T=6C}
\includegraphics[width=0.55\textwidth]{1LineEz_T=10C}\\
\includegraphics[width=0.55\textwidth]{1LineEz_T=20C}
\includegraphics[width=0.55\textwidth]{1LineEz_T=30C}\\
	\caption{ Ez at T = 6, 10, 20 and 30 with $J_x = 10000 \sin(2\pi f_0 t)$ for the left current and $J_x = -10000 \sin(2\pi f_0 t)$ for the right side current.}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{1LineEz_T=10CS}
\includegraphics[width=0.5\textwidth]{1LineEz_T=20CS}\\
\includegraphics[width=0.5\textwidth]{1LineEz_T=30CS}
\includegraphics[width=0.5\textwidth]{1LineEz_T=40CS}\\
	\caption{ Ez at T = 10, 20, 30, 40 with $J_x =  \sin(2\pi f_0 x) \times \sin(2\pi f_0 t)$ for the left current and $J_x = -\sin(2\pi f_0 x)$ for the right side current.}
\end{figure}

Figures 4.11 to 4.12 show at the contour figures of Electric field in the z direction as there is an electric current source moving in the x direction towards the metallic plate from both sides. Figure 4.11 has an electric source current which is only defined in terms of time. As the current starts to move towards the metal and finally meets one another, it starts to move in the y-direction. As time increases, the current moves throughout the center conductor. Figure 4.12 includes an electric current source which moves in both x and t direction and as time increases, the current starts to move towards one another and finally once it reaches, it starts to move in the y-direction. This figure shows similar results to Figure 4.11. The current also moves throughout the center conductor.\\


\begin{figure}	{\textbf{\textit{Simulation 2 Domain}}}
	\begin{center}
		\begin{tikzpicture}
		\filldraw[fill=pink, draw=black] (1,0.3,0)--(1,0.3,4)--(0,0.3,4)--(0,0.3,0)--cycle;
		\filldraw[fill=pink, draw=pink] (0,0.3,4)--(1,0.3,4)--(1,0,4)--(0,0,4)--cycle;
		\filldraw[fill=pink, draw=pink] (1,0.3,0)--(1,0.3,4)--(1,0,4)--(1,0,0)--cycle;
		
		\filldraw[fill=cyan, draw=cyan] (1,0,4)--(3,0,4)--(3,0.3,4)--(1,0.3,4)--cycle;
		\filldraw[fill=cyan, draw=cyan] (1,0.3,4)--(3,0.3,4)--(3,0.3,0)--(1,0.3,0)--cycle;
		
		\filldraw[fill=pink, draw=pink] (3,0.3,0)--(3,0.3,4)--(4,0.3,4)--(4,0.3,0)--cycle;
		\filldraw[fill=pink, draw=pink]  (3,0,4)--(4,0,4)--(4,0.3,4)--(3,0.3,4)--cycle;
		\filldraw[fill=pink, draw=pink] (4,0.3,0)--(4,0.3,4)--(4,0,4)--(4,0,0)--cycle;
		
		\filldraw[fill=cyan, draw=cyan] (0,0.3,0)--(-2,0.3,0)--(-2,0.3,4)--(0,0.3,4)--cycle;
		\filldraw[fill=cyan, draw=cyan] (0,0.3,4)--(-2,0.3,4)--(-2,0,4)--(0,0,4)--cycle;
		
		\filldraw[fill=cyan, draw=cyan] (4,0.3,0)--(6,0.3,0)--(6,0.3,4)--(4,0.3,4)--cycle;
		\filldraw[fill=cyan, draw=cyan] (4,0.3,4)--(6,0.3,4)--(6,0,4)--(4,0,4)--cycle;
		\filldraw[fill=cyan, draw=cyan] (6,0.3,4)--(6,0.3,0)--(6,0,0)--(6,0,4)--cycle;
		
		\draw[thick,->] (6,0,0) -- (6.5,0,0) node[anchor=north east]{\textbf{x}};
		\draw[thick,->] (-2,0.5,0) -- (-2,1,0) node[anchor=north]{\textbf{z}};
		\draw[thick,->] (-2,0,4) -- (-2,0,5) node[anchor=south]{\textbf{y}};
		
		\draw [red,thick,->] (0,0.3,0) -- (1,0.3,0);
		\draw [red,thick,->] (0,0.3,0.5) -- (1,0.3,0.5);
		\draw [red,thick,->] (0,0.3,1) -- (1,0.3,1);
		\draw [red,thick,->] (0,0.3,1.5) -- (1,0.3,1.5);
		\draw [red,thick,->] (0,0.3,2) -- (1,0.3,2);
		\draw [red,thick,->] (0,0.3,2.5) -- (1,0.3,2.5);
		\draw [red,thick,->] (0,0.3,3) -- (1,0.3,3);
		\draw [red,thick,->] (0,0.3,3.5) -- (1,0.3,3.5);
		\draw [red,thick,->] (0,0.3,4) -- (1,0.3,4);
		
		\draw [red,thick,<-] (3,0.3,0) -- (4,0.3,0);
		\draw [red,thick,<-] (3,0.3,0.5) -- (4,0.3,0.5);
		\draw [red,thick,<-] (3,0.3,1) -- (4,0.3,1);
		\draw [red,thick,<-] (3,0.3,1.5) -- (4,0.3,1.5);
		\draw [red,thick,<-] (3,0.3,2) -- (4,0.3,2);
		\draw [red,thick,<-] (3,0.3,2.5) -- (4,0.3,2.5);
		\draw [red,thick,<-] (3,0.3,3) -- (4,0.3,3);
		\draw [red,thick,<-] (3,0.3,3.5) -- (4,0.3,3.5);
		\draw [red,thick,<-] (3,0.3,4) -- (4,0.3,4);
		
		\filldraw [fill={rgb:black,1;white,2}, draw=black, dashed, fill opacity=0.2] (-2,0.3,2) -- (6,0.3,2)--(6,0,2)--(-2,0,2)--cycle;
		\end{tikzpicture} \begin{tikzpicture}
		\draw  (0,0)--(4,0)--(4,1.5)--(0,1.5)--cycle;
		\filldraw[fill=cyan, draw=cyan] (-2,0)--(0,0)--(0,1)--(-2,1)--cycle;
		\filldraw[fill=pink, draw=pink] (0,0)--(1,0)--(1,1)--(0,1)--cycle;
		\filldraw[fill=cyan, draw=cyan] (1,0)--(3,0)--(3,1)--(1,1)--cycle;
		\filldraw[fill=pink, draw=pink] (3,0)--(4,0)--(4,1)--(3,1)--cycle;
		\filldraw[fill=cyan, draw=cyan] (4,0)--(6,0)--(6,1)--(4,1)--cycle;
		
		\draw[step=0.5cm,gray,very thin] (-2,0) grid (6,1.5);
		\draw [red, thick,->] (0,1) -- (1,1);
		\draw [red, thick,<-] (3,1) -- (4,1);
		\draw[thick,->] (6,0) -- (6.5,0) node[anchor=east]{\textbf{x}};
		\draw[thick,->] (-2,1.5) -- (-2,2) node[anchor=north]{\textbf{z}};
		\end{tikzpicture}
	\end{center}
	\caption{Top: Three Metallic plates running parallel to one another with two Dielectric slabs attached in the center along with current running in the x direction, Bottom: 2D Domain of the above. }
\end{figure}	


Figure 4.13 shows the 3D and cross-section domain of a of three parallel metallic plates running with two dielectric slabs in between them. The electric current is added on lumped ports similarly to simulation 1.\\

\begin{figure}
	\includegraphics[width=0.5\textwidth]{3LineEz_T=314C}
	\includegraphics[width=0.5\textwidth]{3LineEz_T=5C}\\
	\includegraphics[width=0.5\textwidth]{3LineEz_T=10C}
	\includegraphics[width=0.5\textwidth]{3LineEz_T=20C}\\
	\includegraphics[width=0.5\textwidth]{3LineEz_T=30C}
	\includegraphics[width=0.5\textwidth]{3LineEz_T=40C}\\
\caption{ Ez at T = 3.14, 5, 10, 20, 30 and 40  with $J_x = 10000 \sin(2\pi f_0 t)$ for the left current and $J_x = -10000 \sin(2\pi f_0 t)$ for the right side current.}
\end{figure}
\begin{figure}
	\includegraphics[width=0.5\textwidth]{3LineEz_T=10CS}
	\includegraphics[width=0.5\textwidth]{3LineEz_T=20CS}\\
	\includegraphics[width=0.5\textwidth]{3LineEz_T=30CS}
	\includegraphics[width=0.5\textwidth]{3LineEz_T=40CS}\\
		\caption{ Ez at T = 10, 20, 30 and 40 with $J_x =  \sin(2\pi f_0 x) \times \sin(2\pi f_0 t)$ for the left current and $J_x = -\sin(2\pi f_0 x)$ for the right side current.}
\end{figure}

Figures 4.14 - 4.15 show the electric field in the z-direction for different times with electric source in terms of time dependent, and both time and space dependent. Figures 4.14-4.15 show similar results to Figures 4.11 -4.12. In time, the electric wave first starts to move in the x direction and once the waves meet, they start to move in the y direction, i.e., towards the center conductor. Also, in each case, as the wave starts to move from dielectric substrate to metallic conductor, the wave starts to change its shape. The electric field in the conductor is stronger than the one on top of a dielectric substrate.\\


\begin{figure}{\textbf{\textit{Simulation 3 Domain}}}
	\begin{center}
		\begin{tikzpicture}
			\filldraw[fill=pink, draw=black] (1,0.3,0)--(1,0.3,4)--(0,0.3,4)--(0,0.3,0)--cycle;
			\filldraw[fill=pink, draw=pink] (0,0.3,4)--(1,0.3,4)--(1,0,4)--(0,0,4)--cycle;
			\filldraw[fill=pink, draw=pink] (1,0.3,0)--(1,0.3,4)--(1,0,4)--(1,0,0)--cycle;
			
			\filldraw[fill=cyan, draw=cyan] (1,0,4)--(3,0,4)--(3,0.3,4)--(1,0.3,4)--cycle;
			\filldraw[fill=cyan, draw=cyan] (1,0.3,4)--(3,0.3,4)--(3,0.3,0)--(1,0.3,0)--cycle;
			
			\filldraw[fill=pink, draw=pink] (3,0.3,0)--(3,0.3,4)--(4,0.3,4)--(4,0.3,0)--cycle;
			\filldraw[fill=pink, draw=pink]  (3,0,4)--(4,0,4)--(4,0.3,4)--(3,0.3,4)--cycle;
			\filldraw[fill=pink, draw=pink] (4,0.3,0)--(4,0.3,4)--(4,0,4)--(4,0,0)--cycle;
			
			\filldraw[fill=cyan, draw=cyan] (0,0.3,0)--(-2,0.3,0)--(-2,0.3,4)--(0,0.3,4)--cycle;
			\filldraw[fill=cyan, draw=cyan] (0,0.3,4)--(-2,0.3,4)--(-2,0,4)--(0,0,4)--cycle;
			
			\filldraw[fill=cyan, draw=cyan] (4,0.3,0)--(6,0.3,0)--(6,0.3,4)--(4,0.3,4)--cycle;
			\filldraw[fill=cyan, draw=cyan] (4,0.3,4)--(6,0.3,4)--(6,0,4)--(4,0,4)--cycle;
			\filldraw[fill=cyan, draw=cyan] (6,0.3,4)--(6,0.3,0)--(6,0,0)--(6,0,4)--cycle;
			
			\filldraw[fill=pink, draw=pink] (-2,0,4)--(6,0,4)--(6,-1,4)--(-2,-1,4)--cycle;
			\filldraw[fill=pink, draw=pink] (6,0,4)--(6,-1,4)--(6,-1,0)--(6,0,0)--cycle;	
		
			\draw[thick,->] (6,-2,0) -- (6.5,-2,0) node[anchor=north east]{\textbf{x}};
			\draw[thick,->] (-2,0.5,0) -- (-2,1,0) node[anchor=north]{\textbf{z}};
			\draw[thick,->] (-2,-1,4) -- (-2,-1,5) node[anchor=south]{\textbf{y}};
			
			\draw [red,thick,->] (0,0.3,0) -- (1,0.3,0);
			\draw [red,thick,->] (0,0.3,0.5) -- (1,0.3,0.5);
			\draw [red,thick,->] (0,0.3,1) -- (1,0.3,1);
			\draw [red,thick,->] (0,0.3,1.5) -- (1,0.3,1.5);
			\draw [red,thick,->] (0,0.3,2) -- (1,0.3,2);
			\draw [red,thick,->] (0,0.3,2.5) -- (1,0.3,2.5);
			\draw [red,thick,->] (0,0.3,3) -- (1,0.3,3);
			\draw [red,thick,->] (0,0.3,3.5) -- (1,0.3,3.5);
			\draw [red,thick,->] (0,0.3,4) -- (1,0.3,4);
			
			\draw [red,thick,<-] (3,0.3,0) -- (4,0.3,0);
			\draw [red,thick,<-] (3,0.3,0.5) -- (4,0.3,0.5);
			\draw [red,thick,<-] (3,0.3,1) -- (4,0.3,1);
			\draw [red,thick,<-] (3,0.3,1.5) -- (4,0.3,1.5);
			\draw [red,thick,<-] (3,0.3,2) -- (4,0.3,2);
			\draw [red,thick,<-] (3,0.3,2.5) -- (4,0.3,2.5);
			\draw [red,thick,<-] (3,0.3,3) -- (4,0.3,3);
			\draw [red,thick,<-] (3,0.3,3.5) -- (4,0.3,3.5);
			\draw [red,thick,<-] (3,0.3,4) -- (4,0.3,4);
		
		\filldraw [fill={rgb:black,1;white,2}, draw=black, dashed, fill opacity=0.2] (-2,0.3,2) -- (6,0.3,2)--(6,-1,2)--(-2,-1,2)--cycle;
		\end{tikzpicture}
		
		\begin{tikzpicture}
		\draw  (0,0)--(4,0)--(4,1.5)--(0,1.5)--cycle;
		\filldraw[fill=cyan, draw=cyan] (-2,0)--(0,0)--(0,1)--(-2,1)--cycle;
		\filldraw[fill=pink, draw=pink] (0,0)--(1,0)--(1,1)--(0,1)--cycle;
		\filldraw[fill=cyan, draw=cyan] (1,0)--(3,0)--(3,1)--(1,1)--cycle;
		\filldraw[fill=pink, draw=pink] (3,0)--(4,0)--(4,1)--(3,1)--cycle;
		\filldraw[fill=cyan, draw=cyan] (4,0)--(6,0)--(6,1)--(4,1)--cycle;
		\filldraw[fill=pink, draw=pink] (-2,0)--(6,0)--(6,-1)--(-2,-1)--cycle;
		
		\draw[step=0.5cm,gray,very thin] (-2,-1) grid (6,1.5);
		\draw [red, thick,->] (0,1) -- (1,1);
		\draw [red, thick,<-] (3,1) -- (4,1);
		\draw[thick,->] (6,-1) -- (6.5,-1) node[anchor=east]{\textbf{x}};
		\draw[thick,->] (-2,2) -- (-2,2.5) node[anchor=north]{\textbf{z}};
		\end{tikzpicture}
		\caption{Top: Three Metallic plates running parallel to one another with two Dielectric slabs attached in the center all on top of a dielectric substrate along with current running in the x direction , Bottom: 2D Domain of the above.}
	\end{center}	
\end{figure}
Figure 4.16 is the simulation domain for three parallel metallic plates running on top of a dielectric slab. Same electric current moving lumped ports are connected to run current towards the metallic plates. Figures 4.17 to 4.18 show similar results to previous simulations. In these figures, it is noticed that as time increases and as electric waves start to move in the y direction, the electric wave starts to disappear as it begins to move into the dielectric slab. Also, as the wave starts to move into a different media, the wave starts to change shape. The simulations are run for a long time, and they confirm that the EM waves die out in dielectric medium.\\

\begin{figure}
	\includegraphics[width=0.5\textwidth]{3LGEz_T=10C}
	\includegraphics[width=0.5\textwidth]{3LGEz_T=20C}\\
	\includegraphics[width=0.5\textwidth]{3LGEz_T=30C}
	\includegraphics[width=0.5\textwidth]{3LGEz_T=40C}\\
	\includegraphics[width=0.5\textwidth]{3LGEz_T=50C}
	\includegraphics[width=0.5\textwidth]{3LGEz_T=60C}\\
	\caption{ Ez at T = 10, 20, 30, 40, 50 and 60 with $J_x = 10000 \sin(2\pi f_0 t)$ for the left current and $J_x = -10000 \sin(2\pi f_0 t)$ for the right side current.}
\end{figure}

\begin{figure}
	\includegraphics[width=0.5\textwidth]{3LGEz_T=10CS}
	\includegraphics[width=0.5\textwidth]{3LGEz_T=20CS}\\
	\includegraphics[width=0.5\textwidth]{3LGEz_T=30CS}
	\includegraphics[width=0.5\textwidth]{3LGEz_T=40CS}\\
	\includegraphics[width=0.5\textwidth]{3LGEz_T=50CS}
	\includegraphics[width=0.5\textwidth]{3LGEz_T=60CS}\\
		\caption{ Ez at T = 10, 20, 30, 40, 50 and 60 with $J_x =  \sin(2\pi f_0 x) \times \sin(2\pi f_0 t)$ for the left current and $J_x = -\sin(2\pi f_0 x)$ for the right side current.}
\end{figure}



\begin{figure} {\textbf{\textit{Simulation 4 Domain}}}
	\begin{center}
		\begin{tikzpicture}
	\filldraw[fill=pink, draw=black] (1,0.3,0)--(1,0.3,4)--(0,0.3,4)--(0,0.3,0)--cycle;
	\filldraw[fill=pink, draw=pink] (0,0.3,4)--(1,0.3,4)--(1,0,4)--(0,0,4)--cycle;
	\filldraw[fill=pink, draw=pink] (1,0.3,0)--(1,0.3,4)--(1,0,4)--(1,0,0)--cycle;
	
	\filldraw[fill=cyan, draw=cyan] (1,0,4)--(3,0,4)--(3,0.3,4)--(1,0.3,4)--cycle;
	\filldraw[fill=cyan, draw=cyan] (1,0.3,4)--(3,0.3,4)--(3,0.3,0)--(1,0.3,0)--cycle;
	
	\filldraw[fill=pink, draw=pink] (3,0.3,0)--(3,0.3,4)--(4,0.3,4)--(4,0.3,0)--cycle;
	\filldraw[fill=pink, draw=pink]  (3,0,4)--(4,0,4)--(4,0.3,4)--(3,0.3,4)--cycle;
	\filldraw[fill=pink, draw=pink] (4,0.3,0)--(4,0.3,4)--(4,0,4)--(4,0,0)--cycle;
	
	\filldraw[fill=cyan, draw=cyan] (0,0.3,0)--(-2,0.3,0)--(-2,0.3,4)--(0,0.3,4)--cycle;
	\filldraw[fill=cyan, draw=cyan] (0,0.3,4)--(-2,0.3,4)--(-2,0,4)--(0,0,4)--cycle;
	
	\filldraw[fill=cyan, draw=cyan] (4,0.3,0)--(6,0.3,0)--(6,0.3,4)--(4,0.3,4)--cycle;
	\filldraw[fill=cyan, draw=cyan] (4,0.3,4)--(6,0.3,4)--(6,0,4)--(4,0,4)--cycle;
	\filldraw[fill=cyan, draw=cyan] (6,0.3,4)--(6,0.3,0)--(6,0,0)--(6,0,4)--cycle;
	
	\filldraw[fill=pink, draw=pink] (-2,0,4)--(6,0,4)--(6,-1,4)--(-2,-1,4)--cycle;
	\filldraw[fill=pink, draw=pink] (6,0,4)--(6,-1,4)--(6,-1,0)--(6,0,0)--cycle;
	
	\filldraw[fill=cyan, draw=cyan] (-2,-1,4)--(6,-1,4)--(6,-1.3,4)--(-2,-1.3,4)--cycle;
	\filldraw[fill=cyan, draw=cyan] (6,-1,4)--(6,-1.3,4)--(6,-1.3,0)--(6,-1,0)--cycle;
	
	
	\draw[thick,->] (6,-2,0) -- (6.5,-2,0) node[anchor=north east]{\textbf{x}};
	\draw[thick,->] (-2,0.5,0) -- (-2,1,0) node[anchor=north]{\textbf{z}};
	\draw[thick,->] (-2,-1,4) -- (-2,-1,5) node[anchor=south]{\textbf{y}};
	
	\draw [red,thick,->] (0,0.3,0) -- (1,0.3,0);
	\draw [red,thick,->] (0,0.3,0.5) -- (1,0.3,0.5);
	\draw [red,thick,->] (0,0.3,1) -- (1,0.3,1);
	\draw [red,thick,->] (0,0.3,1.5) -- (1,0.3,1.5);
	\draw [red,thick,->] (0,0.3,2) -- (1,0.3,2);
	\draw [red,thick,->] (0,0.3,2.5) -- (1,0.3,2.5);
	\draw [red,thick,->] (0,0.3,3) -- (1,0.3,3);
	\draw [red,thick,->] (0,0.3,3.5) -- (1,0.3,3.5);
	\draw [red,thick,->] (0,0.3,4) -- (1,0.3,4);
	
	\draw [red,thick,<-] (3,0.3,0) -- (4,0.3,0);
	\draw [red,thick,<-] (3,0.3,0.5) -- (4,0.3,0.5);
	\draw [red,thick,<-] (3,0.3,1) -- (4,0.3,1);
	\draw [red,thick,<-] (3,0.3,1.5) -- (4,0.3,1.5);
	\draw [red,thick,<-] (3,0.3,2) -- (4,0.3,2);
	\draw [red,thick,<-] (3,0.3,2.5) -- (4,0.3,2.5);
	\draw [red,thick,<-] (3,0.3,3) -- (4,0.3,3);
	\draw [red,thick,<-] (3,0.3,3.5) -- (4,0.3,3.5);
	\draw [red,thick,<-] (3,0.3,4) -- (4,0.3,4);
	
		\filldraw [fill={rgb:black,1;white,2}, draw=black, dashed, fill opacity=0.2] (-2,0.3,2) -- (6,0.3,2)--(6,-1.3,2)--(-2,-1.3,2)--cycle;
		\end{tikzpicture}
		
		\begin{tikzpicture}
		\draw  (0,0)--(4,0)--(4,1.5)--(0,1.5)--cycle;
		\filldraw[fill=cyan, draw=cyan] (-2,0)--(0,0)--(0,1)--(-2,1)--cycle;
		\filldraw[fill=pink, draw=pink] (0,0)--(1,0)--(1,1)--(0,1)--cycle;
		\filldraw[fill=cyan, draw=cyan] (1,0)--(3,0)--(3,1)--(1,1)--cycle;
		\filldraw[fill=pink, draw=pink] (3,0)--(4,0)--(4,1)--(3,1)--cycle;
		\filldraw[fill=cyan, draw=cyan] (4,0)--(6,0)--(6,1)--(4,1)--cycle;
		\filldraw[fill=pink, draw=pink] (-2,0)--(6,0)--(6,-1)--(-2,-1)--cycle;
		\filldraw[fill=cyan, draw=cyan] (-2,-1)--(6,-1)--(6,-2)--(-2,-2)--cycle;
		
		\draw[step=0.5cm,gray,very thin] (-2,-2) grid (6,1.5);
		\draw [red, thick,->] (0,1) -- (1,1);
		\draw [red, thick,<-] (3,1) -- (4,1);
		\draw[thick,->] (6,-2) -- (6.5,-2) node[anchor=east]{\textbf{x}};
		\draw[thick,->] (-2,2) -- (-2,2.5) node[anchor=north]{\textbf{z}};
		\end{tikzpicture}
		\caption{Top: Three Metallic plates running parallel to one another with two Dielectric slabs attached in the center on top a dielectric substrate which is on top of a metallic plate along with current running in the x direction , Bottom: 2D of the above.}
	\end{center}
\end{figure}

 
 
\begin{figure}
	\includegraphics[width=0.5\textwidth]{3LGCEz_T=10C1_}
	\includegraphics[width=0.5\textwidth]{3LGCEz_T=20C1_}\\
	\includegraphics[width=0.5\textwidth]{3LGCEz_T=30C1_}
	\includegraphics[width=0.5\textwidth]{3LGCEz_T=40C1_}\\
	\caption{Ez at T = 10, 20, 30 and 40 with $J_x = 10000 \sin(2\pi f_0 t)$ for the left current and $J_x = -10000 \sin(2\pi f_0 t)$ for the right side current.}
\end{figure}
\begin{figure}
	\includegraphics[width=0.5\textwidth]{3LGCEz_T=10CS}
	\includegraphics[width=0.5\textwidth]{3LGCEz_T=20CS}\\
	\includegraphics[width=0.5\textwidth]{3LGCEz_T=30CS}
	\includegraphics[width=0.5\textwidth]{3LGCEz_T=40CS}\\
	\caption{Ez at T = 10, 20, 30 and 40 with $J_x =  \sin(2\pi f_0 x) \times \sin(2\pi f_0 t)$ for the left current and $J_x = -\sin(2\pi f_0 x)$ for the right side current.}
\end{figure}\\
\\
Figure 4.19 is a 3D and 2D representation of a grounded coplanar waveguide which is mounted on top of a conductor. Figures 4.20 - 4.21 show the electric field in the z-direction for different times for few current sources as mentioned in the above simulations. Since the dielectric substrate slab is thick enough for the electromagnetic waves to die out by the time it reaches the other side of the dielectric substrate, there is no electric wave moving into the conductor at the bottom of the CPW. \\$ $\\

\newpage
\section{Grounded Coplanar Waveguide}
Now, we can continue the same process in a grounded coplanar waveguide. Grounded coplanar waveguide is made of dielectric substrate underneath the coplanar waveguide. It ensures that the electricity produced in the system does not interfere with other components of a printed circuit board. One of common materials used for a substrate is gallium arsenide, GaAs, with a electric permittivity of 12.9. We can use a similar excitation technique on a grounded coplanar waveguide to see how the waves behave inside a grounded coplanar waveguide. Figure 4.22 represents a side view of a grounded coplanar waveguide. \\
		
		
		\begin{figure}
			\begin{center}
			\includegraphics[width=0.5\textwidth]{GCPW}
			\caption{Extended Domain of a Grounded Coplanar Waveguide.}
		\end{center}
		\end{figure}
		
		\begin{figure}
			\begin{center}
				\begin{tikzpicture}			
				\draw[black,thick] (0,0,0) -- (4,0,0)--(4,4,0) --(0,4,0)--(0,0,0)--(0,0,4)--(4,0,4)--(4,0,0)--(0,0,0)--(0,4,0)--(0,4,4)--(0,0,4)--(0,0,0)--(0,4,0)--(4,4,0)--(4,4,4)--(0,4,4)--(0,0,4)--(4,0,4)--(4,4,4);
				
				\filldraw[fill=yellow, draw=black] (1,1.5,0)--(1,2.5,0)--(1,2.5,4)--(1,1.5,4)--cycle;
				\filldraw[fill=yellow, draw=black] (0,2.5,0)--(1,2.5,0)--(1,2.5,4)--(0,2.5,4)--cycle;
				\filldraw[fill=yellow, draw=black] (0,1.5,4)--(1,1.5,4)--(1,2.5,4)--(0,2.5,4)--cycle;
				\node[cylinder, rotate=45+180, fill=pink, draw=black,thick, shape aspect=1, minimum height=2.75cm] at (1.5,2.5,2) {source};
				
				\filldraw[fill=yellow, draw=black] (1.5,1.5,0)--(1.5,2.5,0)--(1.5,2.5,4)--(1.5,1.5,4)--cycle;
				\filldraw[fill=yellow, draw=black] (1.5,2.5,0)--(2.5,2.5,0)--(2.5,2.5,4)--(1.5,2.5,4)--cycle;
				\filldraw[fill=yellow, draw=black] (2.5,2.5,0)--(2.5,1.5,0)--(2.5,1.5,4)--(2.5,2.5,4)--cycle;
				\filldraw[fill=yellow, draw=black] (1.5,2.5,4)--(2.5,2.5,4)--(2.5,1.5,4)--(1.5,1.5,4)--cycle;
				\node[cylinder, rotate=45+180, fill=pink, draw=black,thick, shape aspect=1, minimum height=2.7cm] at (2.75,2.25,1.5) {source};
				
				\filldraw[fill=yellow, draw=black] (3,2.5,0)--(4,2.5,0)--(4,2.5,4)--(3,2.5,4)--cycle;
				\filldraw[fill=yellow, draw=black] (4,2.5,0)--(4,1.5,0)--(4,1.5,4)--(4,2.5,4)--cycle;
				\filldraw[fill=yellow, draw=black] (3,2.5,4)--(4,2.5,4)--(4,1.5,4)--(3,1.5,4)--cycle;
				
				\filldraw [fill=green, draw=black, dashed, fill opacity=0.2] (0,4,2) -- (4,4,2)--(4,0,2)--(0,0,2)--cycle;
				
				
				\draw[thick,->] (5,0,0) -- (6,0,0) node[anchor=north east]{\textbf{x}};
				\draw[thick,->] (0,5,0) -- (0,6,0) node[anchor=north]{\textbf{z}};
				\draw[thick,->] (0,0,5) -- (0,0,6) node[anchor=south]{\textbf{y}};
				
				
				\end{tikzpicture}
				\begin{tikzpicture}
				\draw (0,0)--(4,0)--(4,4)--(0,4)--cycle;
				\draw [fill=yellow] (0,1.5)--(1,1.5)--(1,2.5)--(0,2.5)--cycle;
				\draw [fill=yellow] (1.5,1.5)--(2.5,1.5)--(2.5,2.5)--(1.5,2.5)--cycle;
				\draw [fill=yellow] (3,1.5)--(4,1.5)--(4,2.5)--(3,2.5)--cycle;
				\draw[step=0.5cm,gray,very thin] (0,0) grid (4,4);
				\filldraw[fill=pink, draw=black] (1.25,2) circle (0.15cm);
				\filldraw[fill=pink, draw=black] (2.75,2) circle (0.15cm);
				\filldraw[fill=pink, draw=pink] (-0.5,-0.5)--(0,-0.5)--(0,4.5)--(-0.5,4.5)--cycle;
				\filldraw[fill=pink, draw=pink] (4,-0.5)--(4.5,-0.5)--(4.5,4.5)--(4,4.5)--cycle;
				\filldraw[fill=pink, draw=pink] (0,0)--(0,-0.5)--(4,-0.5)--(4,0)--cycle;
				\filldraw[fill=pink, draw=pink] (0,4)--(0,4.5)--(4,4.5)--(4,4)--cycle;
				\draw[thick,->] (5,0,0) -- (5.5,0,0) node[anchor=north east]{\textbf{x}};
				\draw[thick,->] (0,5,0) -- (0,5.5,0) node[anchor=north]{\textbf{z}};
				\end{tikzpicture}
			\end{center}
			\caption{Grounded Coplanar waveguide with an electric source, $J_y$.} 
				\end{figure}
				In Figure 4.23, the symmetry in the y-direction reduces our system to TM mode with $H_x$, $E_y$ and $H_z$. Figure 4.23 assumes that the domain consists of air and metallic plates. Considering two electric cylindrical shaped sources in the y-direction, the electric source can be defined as
				\begin{equation}
				Jy_1= -\sin(2\pi f_0 t), Jy_2= -\sin(2\pi f_0 t)
				\end{equation}
				at both locations. They are located in the middle of two metallic plates. 

\begin{figure}
	\includegraphics[width=0.5\textwidth]{mesh_Ey_T=1}
	\includegraphics[width=0.5\textwidth]{mesh_Ey_T=2}\\
	\includegraphics[width=0.5\textwidth]{mesh_Ey_T=3}\\
	\caption{Distribution of Ey at T=1, 2, and 3 with sources $Jy_1= -\sin(2\pi f_0 t)$, $Jy_2= -\sin(2\pi f_0 t)$ with air and metallic plates.}
\end{figure}		
Figure 4.24 represents the electric field in the y direction for T=1, 2,and 3. In each of these figures, electric current source is moving in the y-direction. As time increases, the electric field moves towards the metallic plates. However, since the electric permittivitty is lower in air due to the air density, the electric field is stronger in the air domain. \\

\begin{figure}
	\includegraphics[width=0.5\textwidth]{mesh_Ey_s_T=2}
	\includegraphics[width=0.5\textwidth]{mesh_Ey_s_T=3}\\
	\includegraphics[width=0.5\textwidth]{mesh_Ey_s_T=4}
	\includegraphics[width=0.5\textwidth]{mesh_Ey_s_T=42}\\
	\caption{Distribution of $E_y$ at T=2, 3, 4 and 4.2 with sources $Jy_1= -\sin(2\pi f_0 t)$, $Jy_2= -\sin(2\pi f_0 t)$ with dielectric substrate and metallic plates.}
\end{figure}
\begin{figure}
	\includegraphics[width=0.5\textwidth]{contour_Ey_s_T=3}
	\includegraphics[width=0.5\textwidth]{contour_Ey_s_T=4}\\
	\includegraphics[width=0.5\textwidth]{contour_Ey_s_T=42}
	\includegraphics[width=0.5\textwidth]{contour_Ey_s_T=5}\\
	\includegraphics[width=0.5\textwidth]{contour_Ey_s_T=6}\\
	\caption{Contour figures of $E_y$ at T=3, 4, 4.2, 5 and 6 with sources $Jy_1= -\sin(2\pi f_0 t)$, $Jy_2= -\sin(2\pi f_0 t)$ with dielectric substrate and metallic plates.}
\end{figure}

Figure 4.25 shows the mesh of Electric field in the y direction at different times. As time increases, the electric field starts to move into the metallic plates. Eventually, the electric field becomes very strong in the gold plate and vanishes in the dielectric substrate. This is more clear in Figure 4.26. Figure 4.26 presents the contour figures for Electric field at different times. As time increases, the electric field is only present in the metallic plate.


	
\section{Elevated Coplanar Waveguide}
	Paper \cite{Hofschen1996} simulated the results that as the center conductor is elevated by 3$\mu m$, there is strong decrease in the losses of an elevated CPW. The decrease in loss is assumed to be due to more homogeneous field distribution near the slots of the CPW for the elevated center conductor. In the paper, the researcher also proved the phase velocity, $v_p$, is approximately doubled. In this section, we can also look at the behaviour of an electric current in an elevated coplanar waveguide. Figure 4.27 is the schematic view of an elevated CPW from \cite{Hofschen1996}. There is a symmetry in the y-direction, which simplifies the 3D Maxwell's equations to a TE or TM model. For the simulations below, we use a TE model with electric field in the x and z direction and a magnetic field in the y-direction, $E_x$, $E_z$ and $H_y$.
	 \\
	\begin{figure}
	\begin{center}
		\includegraphics[width=0.7\textwidth]{example2}\\
		\caption{Elevated CPW.}
	\end{center}
\end{figure}		
\begin{figure}
	\begin{tikzpicture}
	\draw[step=0.5cm,gray,very thin] (0,-1) grid (5,3);
	\filldraw[fill=cyan, draw=cyan] (0,0)--(1,0)--(1,1)--(0,1)--cycle;
	\filldraw[fill=pink, draw=pink] (1,0)--(4,0)--(4,1)--(1,1)--cycle;
	\filldraw[fill=cyan, draw=cyan] (4,0)--(5,0)--(5,1)--(4,1)--cycle;
	\filldraw[fill=cyan, draw=cyan] (2,1)--(3,1)--(3,2)--(2,2)--cycle;
	\filldraw[fill=pink, draw=pink] (0,0)--(5,0)--(5,-1) --(0,-1)--cycle;
	\filldraw[fill=red, draw=red] (1.3,1)--(1.8,1)--(1.8,1.01)--(1.3,1.01)--cycle;
	\filldraw[fill=red, draw=red] (3.3,1)--(3.8,1)--(3.8,1.01)--(3.3,1.01)--cycle;
\node[draw] at (3,-1.5) {a};

		\draw[step=0.5cm,gray,very thin] (0+8,-1) grid (5+8,3);
		\filldraw[fill=cyan, draw=cyan] (0+8,0)--(1+8,0)--(1+8,1)--(0+8,1)--cycle;
		\filldraw[fill=pink, draw=pink] (1+8,0)--(2+8,0)--(2+8,1)--(1+8,1)--cycle;
			\filldraw[fill=pink, draw=pink] (3+8,0)--(4+8,0)--(4+8,1)--(3+8,1)--cycle;
			\filldraw[fill=pink, draw=pink] (2+8,0)--(3+8,0)--(3+8,1.3)--(2+8,1.3)--cycle;
		\filldraw[fill=cyan, draw=cyan] (4+8,0)--(5+8,0)--(5+8,1)--(4+8,1)--cycle;
		\filldraw[fill=cyan, draw=cyan] (2+8,1+.3)--(3+8,1+.3)--(3+8,2+.3)--(2+8,2+.3)--cycle;
		\filldraw[fill=pink, draw=pink] (0+8,0)--(5+8,0)--(5+8,-1) --(0+8,-1)--cycle;
		\filldraw[fill=red, draw=red] (1.3+8,1)--(1.8+8,1)--(1.8+8,1.01)--(1.3+8,1.01)--cycle;
		\filldraw[fill=red, draw=red] (3.3+8,1)--(3.8+8,1)--(3.8+8,1.01)--(3.3+8,1.01)--cycle;
	\node[draw] at (11,-1.5) {b};
		\end{tikzpicture}
		\caption{Elevated CPW with two different elevation heights.}
\end{figure}

Figure 4.28 shows the 2D domain of an elevated Coplanar Waveguide at two different heights. For each domain, an magnetic current source in the y direction, $J_y$, is placed at top of the dielectric substrates as shown in red.\\ 

\begin{figure}
	\includegraphics[width=0.5\textwidth]{Hy_T=1}
	\includegraphics[width=0.5\textwidth]{Hy_T=2}\\
	\includegraphics[width=0.5\textwidth]{Hy_T=314}
	\caption{Distribution of Magnetic field, $H_y$, with an initial Magnetic source, $J_y$ with an elevated center plate.}
\end{figure}		

\begin{figure}
	\includegraphics[width=0.5\textwidth]{1Hy_T=1}
	\includegraphics[width=0.5\textwidth]{1Hy_T=2}\\
	\includegraphics[width=0.5\textwidth]{1Hy_T=314}
	\caption{Distribution of Magnetic field, $H_y$, with an initial Magnetic source, $J_y$ with an elevated center plate by extra 0.3mm.}
\end{figure}

\begin{figure}
	\begin{center}
\includegraphics[width=0.8\textwidth]{contour_Hy_t=314}
\includegraphics[width=0.8\textwidth]{contour_Hy_t=2}
\caption{Contour figure of $H_y$ at T=3.14 for Elevated CPW (Top) and with extra 0.3mm (Bottom).}
\end{center}
\end{figure}

Figures 4.29 and 4.30 show the magnetic field in the y direction as time goes from 1 to 3.14. As time increases, the magnetic field starts to move, it spreads to all sides. As the two waves meet in the center, the magnetic becomes stronger and starts to move to the top towards the center conductor. As the wave moves to the side conductor, the wave starts to change shape and the magnetic field again start to become stronger. Figure 4.30, T=3.14 shows more clearly where the magnetic field starts to change as it moves to the metallic plate. As time increases, the magnetic field does continue to move towards the center more elevated conductor. Figure 4.31 show the contour figures of magnetic field in the y-direction at T=3.14 for both cases of elevated height. The top picture in Figure 4.31 shows the magnetic field with less elevated height. In this case, the center conductor is not as high, hence there is a stronger magnetic field around the mesh points (175,225). Comparing it with the bottom picture, the magnetic field is not as strong in the center since the center conductor is elevated slightly more.
\newpage
\singlespacing
\section{Transition between a Coplanar Waveguide and a Rectangular Waveguide}
\doublespacing
Now, we can finally apply the EO-4th-EC scheme to the transition problem in \cite{Deslandes2001}. The structure includes a rectangular waveguide synthesized with a linear array of metalized via holes on the same substrate used for the planar circuit as shown in Figure 4.3. Such integrating structure reduces the size and cost. This structure is made of two back to back transitions from CPW to rectangular waveguide mounted on top same substrate layer. It consists of a coplanar waveguide with 90 degrees bend on each slot. A stub is added on the CPW lines to match the transition and the rectangular waveguide. However, the domain has been slightly modified to as seen in Figure 4.34 which has been even more simplified as seen in Figure 4.32. The same dimensions were used, however, less via metallic holes were considered to simplify the complex structure. In the experimental results concurred by paper \cite{Deslandes2001}, CPW was operated at a frequency of 19 to 38 GHz and 27.5 to 29.5 GHz provided the optimal result in terms of return loss. For this example, 28GHz is considered. The transition has a symmetry about the z-axis, hence, TE mode with $E_x$, $E_y$ and $H_z$ is used. The magnetic current source is added in the center of the domain as seen as red circle in Figure 4.32.\\


\begin{figure}
	\begin{center}
		\begin{turn}{90}
		\begin{tikzpicture}
		\draw[step=0.5cm,very thin] (1-1,3-3) grid (8-1,8-3);
		\filldraw[fill=cyan, draw=black] (1-1,6.3-3)--(3-1,6.3-3)--(3-1,7.75-3)--(3.3-1,7.75-3)--(3.3-1,5.8-3)--(2.5-1,5.8-3)--(2.5-1,5.9-3)--(1-1,5.9-3)--cycle;
		\filldraw[fill=cyan, draw=black] (1-1,5-3)--(2.5-1,5-3)--(2.5-1,5.1-3)--(3.3-1,5.1-3)--(3.3-1,3.25-3)--(3-1,3.25-3)--(3-1,4.6-3)--(1-1,4.6-3)--cycle;
		\filldraw[fill=cyan, draw=black]
		(8-1,5-3)--(6.5-1,5-3)--(6.5-1,5.1-3)--(5.7-1,5.1-3)--(5.7-1,3.25-3)--(6-1,3.25-3)--(6-1,4.6-3)--(8-1,4.6-3)--cycle;
		\filldraw[fill=cyan, draw=black]
		(8-1,5.9-3)--(6.5-1,5.9-3)--(6.5-1,5.8-3)--(5.7-1,5.8-3)--(5.7-1,7.75-3)--(6-1,7.75-3)--(6-1,6.3-3)--(8-1,6.3-3)--cycle;
			\draw[red,fill=red] (4.5-1,5.5-3) circle (0.3);
		\end{tikzpicture}
		\end{turn}
	\end{center}
	\caption{Cross section of transition between CPW and Rectangular Waveguide with source located in the center of the domain.}
	\end{figure}


\begin{figure}

\includegraphics[width=0.5\textwidth]{HzatT=5_Legs}
\includegraphics[width=0.5\textwidth]{HzatT=10_Legs}\\
\includegraphics[width=0.5\textwidth]{HzatT=15_Legs}
\includegraphics[width=0.5\textwidth]{HzatT=20_Legs}\\
\includegraphics[width=0.5\textwidth]{HzatT=25_Legs}
\includegraphics[width=0.5\textwidth]{HzatT=30_Legs}\\
	\caption{Propagations of $H_z$ at T= 5, 10, 15, 20, 25 and 30. }
\end{figure}
\newpage
Figure 4.33 shows the contour figures of magnetic field in the z-direction for T=5, 10, 15, 20, 25 and 30. As time increases, the wave starts to move outwards towards the metallic stubs. Eventually the magnetic field is stronger at the four corners of the metal and starts to move in the four legs of the transition.\\ 

\begin{figure}
	\begin{center}
		\begin{turn}{90}
		\begin{tikzpicture}
		\draw[step=0.5cm,very thin] (1-1,3-3) grid (8-1,8-3);
		\filldraw[fill=cyan, draw=black] (1-1,6.3-3)--(3-1,6.3-3)--(3-1,7.75-3)--(3.3-1,7.75-3)--(3.3-1,5.8-3)--(2.5-1,5.8-3)--(2.5-1,5.9-3)--(1-1,5.9-3)--cycle;
		\filldraw[fill=cyan, draw=black] (1-1,5-3)--(2.5-1,5-3)--(2.5-1,5.1-3)--(3.3-1,5.1-3)--(3.3-1,3.25-3)--(3-1,3.25-3)--(3-1,4.6-3)--(1-1,4.6-3)--cycle;
		\filldraw[fill=cyan, draw=black]
		(8-1,5-3)--(6.5-1,5-3)--(6.5-1,5.1-3)--(5.7-1,5.1-3)--(5.7-1,3.25-3)--(6-1,3.25-3)--(6-1,4.6-3)--(8-1,4.6-3)--cycle;
		\filldraw[fill=cyan, draw=black]
		(8-1,5.9-3)--(6.5-1,5.9-3)--(6.5-1,5.8-3)--(5.7-1,5.8-3)--(5.7-1,7.75-3)--(6-1,7.75-3)--(6-1,6.3-3)--(8-1,6.3-3)--cycle;
		\draw[black,fill=cyan] (2.5-1,7-3) circle (0.2);
		\draw[black,fill=cyan] (2.5-1,4-3) circle (0.2);
		\draw[black,fill=cyan] (2.5+2-1,7-3) circle (0.2);
		\draw[black,fill=cyan] (2.5+2-1,4-3) circle (0.2);
		\draw[black,fill=cyan] (2.5+4-1,7-3) circle (0.2);
		\draw[black,fill=cyan] (2.5+4-1,4-3) circle (0.2);
		\draw[red,fill=red] (4.5-1,5.5-3) circle (0.3);
		\end{tikzpicture}
	\end{turn}
	\end{center}
	\caption{Cross section of transition between CPW and Rectangular Waveguide with source located in the center of the domain.}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{HzatT=5C_LegsCircle}
\includegraphics[width=0.5\textwidth]{HzatT=10C_LegsCircle}\\
\includegraphics[width=0.5\textwidth]{HzatT=15C_LegsCircle}
\includegraphics[width=0.5\textwidth]{HzatT=20C_LegsCircle}\\
\includegraphics[width=0.5\textwidth]{HzatT=25C_LegsCircle}
\includegraphics[width=0.5\textwidth]{HzatT=27C_LegsCircle}\\

	\caption{Propagations of $H_z$ at T= 5, 10, 15, 20, 25 and 27 in the presence of via holes.}
\end{figure}


Figure 4.34 shows the cross-section of the transition between CPW and a rectangular waveguide in the presence of the via holes. Figure 4.35 shows the contour figures of magnetic field in the z-direction for T = 5, 10, 15, 20, 25 and 27. As time increases, the wave starts to move outwards towards the metallic stubs and the via holes. The dielectric substrate has a higher relative permittivitty than the material used in the metals and holes. Hence, the electric field is stronger where there is a lower electric permittivitty. The magnetic field is symmetric throughout the grid. Via holes are located near mesh points (125,125),(125,175), (125,275) (275, 125), (275,175) and (275,275). In Figure 4.35, the bottom left picture show the magnetic field at the location of the holes at mesh point (125,175) and (275,175). When T=10, the magnetic field hits the four corners and the wave starts to propagate differently. As the wave hits the holes, those sections are blue colour as seen in the last picture of Figure 4.35.\\
 


\begin{figure}
	\begin{center}
		\begin{turn}{90}
		\begin{tikzpicture}
		\draw[step=0.5cm,very thin] (1-1,3-3) grid (8-1,8-3);
		\filldraw[fill=cyan, draw=black] (1-1,6.3-3)--(3-1,6.3-3)--(3-1,7.75-3)--(3.3-1,7.75-3)--(3.3-1,5.8-3)--(2.5-1,5.8-3)--(2.5-1,5.9-3)--(1-1,5.9-3)--cycle;
		\filldraw[fill=cyan, draw=black] (1-1,5-3)--(2.5-1,5-3)--(2.5-1,5.1-3)--(3.3-1,5.1-3)--(3.3-1,3.25-3)--(3-1,3.25-3)--(3-1,4.6-3)--(1-1,4.6-3)--cycle;
		\filldraw[fill=cyan, draw=black]
		(8-1,5-3)--(6.5-1,5-3)--(6.5-1,5.1-3)--(5.7-1,5.1-3)--(5.7-1,3.25-3)--(6-1,3.25-3)--(6-1,4.6-3)--(8-1,4.6-3)--cycle;
		\filldraw[fill=cyan, draw=black]
		(8-1,5.9-3)--(6.5-1,5.9-3)--(6.5-1,5.8-3)--(5.7-1,5.8-3)--(5.7-1,7.75-3)--(6-1,7.75-3)--(6-1,6.3-3)--(8-1,6.3-3)--cycle;
		\draw[black,fill=cyan] (2.5-1,7-3) circle (0.2);
		\draw[black,fill=cyan] (2.5-1,4-3) circle (0.2);
		\draw[black,fill=cyan] (2.5+2-1,7-3) circle (0.2);
		\draw[black,fill=cyan] (2.5+2-1,4-3) circle (0.2);
		\draw[black,fill=cyan] (2.5+4-1,7-3) circle (0.2);
		\draw[black,fill=cyan] (2.5+4-1,4-3) circle (0.2);
		\draw[red,fill=red] (2,5.5-3) circle (0.2);
		\draw[red,fill=red] (5,5.5-3) circle (0.2);
		\end{tikzpicture}
		\end{turn}
	\end{center}
	\caption{Cross section of transition between CPW and Rectangular Waveguide with two sources.}
\end{figure}
\newpage
Figure 4.36 shows the cross-section of the transition similar to Figure 4.34, but this transition is excited at two different points with magnetic current sources. Figure 4.37 shows the contour figures of the magnetic field in the z-direction for different times when there are two sources present. In the first picture, the magnetic field, $H_z$, starts to approach one another. At T=10, the magnetic field lines intersect one another and by T=15, the magnetic field has reached the four metallic corners. At this time, wave has just reached the holes located at mesh points (125,175) and (275,175). It is clear that at T = 20, the wave changes its shape as it makes two almost perfect circles at via holes. The wave continues to pass and move through the media with an increase in time.

\begin{figure}

	\includegraphics[width=0.5\textwidth]{2SHzatT=5}
	\includegraphics[width=0.5\textwidth]{2SHzatT=10}\\
	\includegraphics[width=0.5\textwidth]{2SHzatT=15}
	\includegraphics[width=0.5\textwidth]{2SHzatT=20}\\
	\includegraphics[width=0.5\textwidth]{2SHzatT=25}
	\includegraphics[width=0.5\textwidth]{2SHzatT=27}\\
	
	\caption{Propagation of $H_z$ at T= 5, 10, 15, 20, 25 and 27 in the presence of via holes and two sources.}
\end{figure}


\newpage

\chapter{Conclusion}

In this thesis, we developed a new even-odd cycled energy-conserved splitting FDTD scheme to solve Maxwell's equations in 2D. The scheme is second order in time and fourth order in space. Theoretical analysis was carried out to confirm that the scheme is unconditionally stable and non-dissipative.
Further, more theoretical analysis was taken to prove the convergence of the scheme and to show that the scheme is energy-conserved. Numerical experiments were computed using MATLAB codes to confirm that the scheme conserves energy, and has high order accuracy in space and time. \\
  
The developed EO-4th-EC method can be applied to solve high frequency problems in many applications. The method was applied to compute the electromagnetic field in coplanar waveguides. Coplanar waveguides are usually a part of a printed circuit boards found in many integrated circuits. Two lumped ports are added to excite a coplanar waveguide where the current is either moving towards or away from the center conductor in each port. Further, a grounded coplanar waveguide  and an elevated coplanar waveguide was modeled. Last experiment of Chapter 4 looked at a transition between a coplanar waveguide and a rectangular waveguide. This transition is part of a printed circuit board which is proved to be more optimal in terms of size and cost. Although, Gaussian pulse was initially considered, Sinusoidal sources were used throughout to model each scenario. \\

As for the future work, the EO-4th-EC scheme can be modified and be used for all sort of transmission lines and can be used to compute the EM waves in coaxial cable, microstrip, striplines or any transitions between them to name a few. The proposed scheme can further be extended to 3D Maxwell's equations. However, more research, analysis and simulations are needed. 





\newpage
\begin{thebibliography}{9}
	
\bibitem{Adve1997}	
Adve, R., \textit{``Extrapolation of time-domain responses from three-dimensional conducting objects utilizing the matrix pencil technique"}. IEEE Transactions on Antennas and Propagation 45.1 (1997): 147-156.	
	
\bibitem{Benavides2011}	
Benavides, M., \emph{``A novel set of reduced equations to model perfect layer matched (PML) in FDTD."} Revista mexicana de física E 57.1 (2011): 25-31.
	
\bibitem{Berenger1994}
Berenger, J.P., \emph{``A perfectly matched layer for the absorption of electromagnetic waves."} Journal of computational physics 114.2 (1994): 185-200.

\bibitem{Cangellaris1993}
Cangellaris, A. C., \textit{``Numerical stability and numerical dispersion of a compact 2-D/FDTD method used for the dispersion analysis of waveguides.}" IEEE Microwave and Guided Wave Letters 3.1 (1993): 3-5.

\bibitem{Chen2008}
Chen, W., Li, X., and Liang, D., \emph{``Energy-conserved splitting FDTD methods for Maxwell's equations"}, Numerische Mathematik, 108.3 (2008): 445-485.

\bibitem{Deslandes2001}
Deslandes, D., Wu, K., \emph{``Integrated Transition of Coplanar to Rectangular Waveguides"}, Microwave Symposium Digest, 2001 IEEE MTT-S International, Vol. 2. IEEE (2001): 619-622.

\bibitem{Gaffar2014}
Gaffar, M., and Dan J., \textit{``An explicit and unconditionally stable FDTD method for electromagnetic analysis}." IEEE Transactions on Microwave Theory and Techniques 62.11 (2014): 2538-2550.

\bibitem{gedney1996}
Gedney, S. D., \emph{``An anisotropic PML absorbing media for the FDTD simulation of fields in lossy and dispersive media."} Electromagnetics 16.4 (1996): 399-415.

\bibitem{Guo2004}
Guo, S., \textit{``Loss and dispersion analysis of microstructured fibers by finite-difference method.}" Optics Express 12.15 (2004): 3341-3352.

\bibitem{Hirono2000}
Hirono, T., \emph{``The second-order condition for the dielectric interface orthogonal to the Yee-lattice axis in the FDTD scheme."} IEEE microwave and guided wave letters 10.9 (2000): 359-361.

\bibitem{Hofschen1996}
Hofschen, S., and Wolff, I., \emph{``Simulation of an elevated coplanar waveguide using 2-D FDTD}", Microwave and Guided Wave Letters, IEEE, 6.1 (1996): 28-30.

\bibitem{Hong2014}
Hong, J., Lihai J., and Kong, L., \emph{``Energy-dissipation splitting finite-difference time-domain method for Maxwell equations with perfectly matched layers."} Journal of Computational Physics 269 (2014): 201-214.

\bibitem{hwang2001}
Hwang, K., and Andreas C. C., \emph{``Effective permittivities for second-order accurate FDTD equations at dielectric interfaces."} IEEE Microwave and Wireless Components Letters 11.4 (2001): 158-160.

\bibitem{Inan2011}
Inan, U. and Marshall, R., \emph{``Numerical Electromagnetics; The FDTD Method}", Cambridge University, New York, 2011.

\bibitem{Jackson1986}
Jackson, R. W., \emph{``Considerations in the use of coplanar waveguide for millimeter-wave integrated circuits."} IEEE Transactions on microwave theory and techniques 34.12 (1986): 1450-1456.

\bibitem{Kunz1993}
Kunz, K. S., and Luebbers, R. J., \emph{The finite difference time domain method for electromagnetics.} CRC press, 1993.

\bibitem{Lee2003}
Lee, J., and Fornberg, B., ``A split step approach for the 3-D Maxwell's equations." Journal of Computational and Applied Mathematics 158.2 (2003): 485-505.

\bibitem{lee2008}
Lee, S., Jung, S., and Lee, H., \emph{``Ultra-wideband CPW-to-substrate integrated waveguide transition using an elevated-CPW section."} IEEE Microwave and Wireless Components Letters 18.11 (2008): 746-748.

\bibitem{Li2015}
Li, W. S., and Liang, D., \emph{``Modeling dielectric waveguides and PML scatterings by the efficient EC-S-FDTD scheme."} Control, Mechatronics and Automation Technology: Proceedings of the International Conference on Control, Mechatronics and Automation Technology (ICCMAT 2014), July 24-25, 2014, Beijing, China. Vol. 6. CRC Press, 2015: 227-230.

\bibitem{Li2016}
Li, W., and Liang, D., \textit{``Symmetric Energy-Conserved S-FDTD Scheme for Two-Dimensional Maxwell’s Equations in Negative Index Metamaterials}." Journal of Scientific Computing 69.2 (2016): 696-735.


\bibitem{Liang2013}
Liang, D., and Yuan, Q., \emph{``The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations}", Journal of Computational Physics, 243 (2013): 344-364.

\bibitem{Liang1989}
Liang, G-C., Liu, Y-U., and Mei,K.K., \emph{``Full-wave analysis of coplanar waveguide and slotline using the time-domain finite-difference method"}, Microwave Theory and Techniques, IEEE Transactions  37.12 (1989): 1949-1957.

\bibitem{Loui2004}
Loui, H., \emph{``1D-FDTD using MATLAB."} ECEN-6006 Numerical Methods in Photonics Project-1 (2004): 1-13.

\bibitem{Marchuk1990}
Marchuk, G-I., \emph{``Splitting and alternating direction methods}", Handbook of numerical analysis. 1 (1990).

\bibitem{Maloney1992}
Maloney, J. G., and Smith, G.S., \emph{``The efficient modeling of thin material sheets in the finite-difference time-domain (FDTD) method."} IEEE Transactions on antennas and Propagation 40.3 (1992): 323-330.

\bibitem{Monk1994}
Monk, P. and Suli E., \emph{``A convergence analysis of Yee's scheme on nonuniform grids}", SIAM Journal of Numerical Analysis, 31 (1994), 1056-1070.

\bibitem{Namiki1999}
Namiki, T., \emph{``A new FDTD algorithm based on alternating-direction implicit method}", Microwave Theory and Techniques, IEEE Transactions, 47.10 (1999): 2003-2007.

\bibitem{Nguyen2015}
Nguyen, B T., Furse, C., and Simpson, J. J., \textit{``A 3-D stochastic FDTD model of electromagnetic wave propagation in magnetized ionosphere plasma}." IEEE Transactions on Antennas and Propagation 63.1 (2015): 304-313.

\bibitem{ponchak1990}
Ponchak, G. E., and Simons, R.N., \emph{``A new rectangular waveguide to coplanar waveguide transition."} (1990).

\bibitem{Schneider2010}
Schneider, J. B., \emph{``Understanding the finite-difference time-domain method."} School of electrical engineering and computer science Washington State University. (2010).

\bibitem{simons2004}
Simons, R. N., \emph{Coplanar waveguide circuits, components, and systems.} Vol. 165. John Wiley $\&$ Sons, 2004.

\bibitem{Singh2011}
Singh, P., and Verma, A. K., \emph{``Dielectric loss computation of multilayer Coplanar Waveguide."} Communications and Signal Processing (ICCSP), 2011 International Conference on. IEEE, 2011:428-430.

\bibitem{Sullivan2013}
Sullivan, D. M., \emph{Electromagnetic simulation using the FDTD method.} John Wiley $\&$ Sons, 2013.

\bibitem{Taflove2000}
Taflove, A. and Hagness, S., \emph{``Computational Electromagnetics; The Finite Difference Time-Domain Method}", Second ed., Artech House, Boston, MA, 2000.

\bibitem{toroglu2014}
Toroğlu, G., and Sevgi, L., \emph{``Finite-difference time-domain (FDTD) MATLAB codes for first-and second-order em differential equations."} (2014):221-239.

\bibitem{verma2005}
Verma, A. K., and Singh, H., \emph{``Dielectric loss of multilayer coplanar waveguide using the single layer reduction (SLR) formulation."} 2005 Asia-Pacific Microwave Conference Proceedings. Vol. 1. IEEE, 2005: 3-pp.

\bibitem{verma2012}
Verma, A. K., and Singh, P., \emph{``Loss Computation of Multilayer Coplanar Waveguide using Single Layer Reduction Method."} International Journal of Electromagnetics and Applications 2.6 (2012): 174-181.

\bibitem{welfer2007}
Welfert, B.D., `\textit{`Analysis of Iterated ADI-FDTD schemes for Maxwell's equations.}" J. Comp. Phys., 222(2007),9-27.

\bibitem{Xie2002}
Xie, Z., Chan, C-H., and Zhang, B., \emph{``An explicit fourth-order staggered finite-difference time-domain method for Maxwell's equations}", Journal of Computational and Applied Mathematics, 147.1 (2002): 75-98.

\bibitem{xu2005}
Xu, F., and Wu, K., \emph{``Guided-wave and leakage characteristics of substrate integrated waveguide."} IEEE Transactions on microwave theory and techniques 53.1 (2005): 66-73.

\bibitem{Yanenko1971}
Yanenko, N. N., \emph{``The method of fractional steps}". Springer-Verlag, New York, 1971.

\bibitem{Yee1966}
Yee, K. S., \emph{``Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media}", IEEE Trans. Antennas Propag, 14.3 (1966): 302-307.

\bibitem{Yefet2001}
Yefet, A., and Petropoulos, P. G., \textit{ ``A staggered fourth-order accurate explicit finite difference scheme for the time-domain Maxwell's equations."} Journal of Computational Physics 168.2 (2001): 286-315.

\bibitem{Zhao2002}
Zhao, A. P., \emph{``Analysis of the numerical dispersion of the 2D alternating-direction implicit FDTD method}", Microwave Theory and Techniques, IEEE Transactions, 50.4 (2002): 1156-1164.

\bibitem{Zhao2016}
Zhao, L., Li, W., and Yu, W., \textit{ ``Numerical dispersion analysis of energy conserved splitting FDTD method for Maxwell's equations."} Computational Electromagnetics (ICCEM), 2016 IEEE International Conference on. IEEE, 2016: 333-335.

\bibitem{Zheng2000}
Zheng, F., Chen, Z. and Zhang, J., \emph{``Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method}", IEEE Trans, Microwave Theory Tech., 48(2000): 1550-1558.
\end{thebibliography}
\end{document}
