# Photonic Project: FDTD Simulation of Waveguide Splitter
contributed by < WeiTing Chien [tintinjian12999](https://github.com/tintinjian12999) >
#### For .gif file to work well, view this file on https://hackmd.io/@tintinjian12999/PhotonicProject
## Introduction
Owing to the exponential-growing data-rate requirement. The next generation communication protocol (6G) aiming the communication frequency band to THz (0.1THz ~ 10THz). Which has wide unallocated bands and is essential for high datarate (Tbit/s) transmission.

Traditional electrical approaches on communications system, with adventages on high signal power and quality, compatible size, and mature manifacture technique on integrated circuits (ICs). However, they suffers from high attenuation loss at THz frequency region and the operation bandwidth limitation on electrical devices.

On the the other hand, the photonic approaches, with much wider operation bandwidths, and much smaller power attentuation. Has became a hit.

Moreover, with the progress on silicon photonics which integrated with the CMOS process, it is possible to fabricate low-cost, tiny-scale, and mass-produce the photonic integrated circuits(PICs), which transmit signals through lights on optical waveguides rather than electrical signals on wires.

Here, in this report, we will studying the fundemental theorys of optical waveguides, and simulate some waveguide splitters devices through Finite-difference time-domain (FDTD) method with the help of a powerful python interface [meep](https://meep.readthedocs.io/en/master/).
## Photonic waveguide

### Waveguide modes
Consider a planar dielectric waveguide with infinite length along x-axis, waveguide width $d$ along y-axis as shown in figure below


The waveguide is consist of two regions: The guided region with refractive index $n1$, and the unguided region with refractive index $n2$, where n1 > n2 in order to fullfill the total internal reflection condition.

To maintain the propagation in waaveguide, the self-consistence condition must be satisfied, As shown in the graph below, the phase shift between the origional wave and the double reflected wave should be multiple of $2 \pi$
leads to

$$
\frac{2\pi}{\lambda}2dsin\theta - 2 \phi_r = 2 \pi m
$$

where $\phi_r$ represent the phase shift introduced at the dielectric boundary under the total internal condition.
Under TE condition:
$$
tan \frac{\phi_r}{2} = \sqrt{\frac{cos^2\theta_c}{cos^2\bar{\theta}}-1}
$$
TM:
$$
tan \frac{\phi_r}{2} = \frac{-1}{sin^2\theta_c} \sqrt{\frac{cos^2\theta_c}{cos^2\bar{\theta}}-1}
$$

where $\theta_c$ represents the critical angle.

Consider TE case, substitute the waveguide phase shift condition of the waveguide into the equation above, it than be:

$$
tan(\pi \frac{d}{\lambda}sin\theta - m \frac{\pi}{2}) = \sqrt{\frac{sin^2 \bar{\theta}_c}{sin^2\theta} - 1}
$$

Solving by graphic, the number of Mode can than be obtained, as the garphic shown below.
In this case, the number of mode is 9, where the intersection correpoend to $sin\theta_m$ of each propagating modes.

The propagation constant along z direction is
$$
k_z = \beta_m = n_1k_0cos\theta_m
$$


Finally, substituting $cos \bar{\theta}_c = \frac{n_2}{n_1}$, the mode number $M$ can than be described as

$$
M \doteq \frac{2d}{\lambda_0}NA
$$
where $NA = \sqrt{n^2_1 - n^2_2}$ and $\doteq$ means increase to the nearest integer.

The mode cutoff frequency which represents the highest frequency for single-mode condition to operate is
$$
v_c = \frac{\omega_c}{2\pi} = \frac{1}{NA}\frac{c_0}{2d}
$$
### Field distribution
Consider Maxwell's equations in medium
$$
\nabla \times \bar{\textbf{E}} = - \mu_0 \frac{\partial \bar{\textbf{H}}}{\partial t}
$$

$$
\nabla \times \bar{\textbf{H}} = \epsilon_0 n^2 \frac{\partial \bar{\textbf{E}}}{\partial t}
$$
Assume the fields take the form
$$
\bar{\textbf{E}} = \textbf{E}(x, y)\exp(j(\omega t - \beta z))
$$
In slab waveguide, the derivative of fields in x direction is zero. The Maxwell equations can be simplified as:

$$
\begin{cases}
\frac{\partial E_z}{\partial y} + j \beta E_y = -j\omega\mu_0H_x \\
-j \beta E_x = -j\omega\mu_0H_y\\
-\frac{\partial E_x}{\partial y} = -j\omega\mu_0H_z
\end{cases}
$$

$$
\begin{cases}
\frac{\partial H_z}{\partial y} + j \beta H_y = -j\omega\epsilon_0n^2E_x \\
-j \beta H_x = -j\omega\epsilon_0n^2E_y\\
-\frac{\partial H_x}{\partial y} = -j\omega\epsilon_0n^2E_z
\end{cases}
$$

Further consider the TE condition, where $E_z = 0$

$$
H_y = \frac{\beta}{\omega\mu_0} E_x
$$

$$
H_z = \frac{1}{j\omega\mu_0} \frac{\partial E_x}{\partial y}
$$

$$
\frac{1}{j\omega\mu_0} \frac{\partial^2 E_x}{\partial^2 y} + j\beta^2 \frac{1}{\omega\mu_0} E_x = -j\omega\epsilon_0n^2E_x
$$

Leads to the final result, the wave equations of $E_x$ for the TE modes:
$$
\frac{d^2E_x}{dy^2} + (k^2n^2 - \beta^2)E_x = 0
$$

let $\beta = \beta_m$ represents propagation constants in defferent modes and $\gamma_m = k^2n^2 - \beta^2$.

In the guided region where $n = n_1$, and $k^2n^2_1 > \beta^2$, $\gamma_m > 0$. It leads to the solution form of the wave equation to be:

$$
E_x(y) = C_1cos(\gamma_my) + C_2sin(\gamma_my)
$$

It represents the fields in the guided region is indeed sinusoidal wave and does not vanish (In lossless case).

On the other hand, in the non-guided region where $n = n_2$, $k^2n^2_2 < \beta^2$, $\gamma_m < 0$, the solution form of the wave equation become:

$$
E_x(y) = C_1\exp(\gamma_my) + C_2exp(-\gamma_my)
$$

Which will decade along the y-axis and represents the so-called evanescent wave.

Similar result in the TM case can be devloped through the steps mentioned above.
## Finite-Difference Time-Domain method (FDTD)
### Introduction
Finite-Difference Time-Domain method or so-called FDTD method is a numerical method aiming to solve differential equations problems. And is widely used to describe electrodynamics topics. It can shows how the field propagate in the space throughout time.
### Discretization

The basic idea of FDTD is the discretization both on space and time. For the space part, discretization is done by the well-known Yee's space lattice:

As illsturated above, the electric fields and magnetic fields are located on difference grids in the space. Also each components will never overlap with each other. Making the numerical more robust and can describing all the components simultaneously.

Meanwhile, on the time part, to fulfill Farady's law and Ampere's law. The electric fields are first updated on time steps. And based on the previously-computed electric fields, the magnetic fields are computed on half-time steps. The magnetic fields is than used to update the electric fields. As the figure shows below.

The origional continuous diffraction equation is discretized to:
$$
\frac{\partial u}{\partial x}(i \Delta x, j \Delta y, k \Delta z, n \Delta t) = \frac{u^n_{i + 1,j,k}-u^n_{i,j,k}}{\Delta x} + O[(\Delta x)^2]
$$

$$
\frac{\partial u}{\partial t}(i \Delta x, j \Delta y, k \Delta z, n \Delta t) = \frac{u^{n+\frac{1}{2}}_{i ,j,k}-u^{n-\frac{1}{2}}_{i ,j,k}}{\Delta t} + O[(\Delta t)^2]
$$

For example, the origional Maxwell equation
$$
\frac{\partial E_x}{\partial t} = \frac{1}{\epsilon}[\frac{\partial H_x}{\partial y} - \frac{\partial H_y}{\partial z}]
$$

Can be discretized as
$$
\frac{E_x|^{n+1/2}_{i, j + 1/2, k + 1/2} - E_x|^{n-1/2}_{i, j + 1/2, k + 1/2}}{\Delta t} = \frac{1}{\epsilon_{i,j+1/2,k+1/2}}[
\frac{H_z|^{n}_{i, j + 1, k + 1/2} - H_z|^{n}_{i, j, k + 1/2}}{\Delta y} - \frac{H_y|^{n}_{i, j + 1/2, k + 1} - H_y|^{n}_{i, j + 1/2, k }}{\Delta z} ]
$$

Obviously, to implement FDTD by our own might be time-consuming and with certain difficulty. Also it is not our main task here. Here, by the help of [meep](https://meep.readthedocs.io/en/master/), a FDTD python development interface, we can customized our own FDTD simulation much more easier.

## Waveguide splitter
A 50-50 wave guide splitter is a common waveguide element that can be used to divided signal from one input to multiple outout. Here we simulate several design on waveguide splitter and see the effects.
### Simulation set up
The wavelength of light is set to be $1.55\mu m$. The waveguide material is assumed to be lossless Si covered with air. The width of wave guide is adjust such that single mode condition is acheived under $\lambda = 1.55\mu m$. The light source is a monochromatic continuous wave which will keep oscilate for 50 time units. The total simulation time will be 400 time units. The Poynting factor sensors which is used to calculate the power of electric field flow through certain area will be used to determine the effect of the splitters. All code can be found in my github [Photonic_Project](https://github.com/tintinjian12999/Photonic_project/tree/main).
### Vertical Splitter
We first observe a intuitive example
The light is simply hit a perpendicular wall and divided into two pass. The blue line represents the Poynting factor sensors.

The resulting $power/m^2$ is

1.187 1.183 0.141 0.141 0.109 0.109 0.095 0.095

The waveguide did propagte well, also from the result of sensor5(S5), S6, S7, S8. The splitter did quite a good work.
However, by observing S2, S3, S4 we can find a huge loss of
$$
\frac{1.183 - 0.141 - 0.141}{1.183} = 76.16\%
$$
Also by observing S3, S5, S4, S6. There are also losses of

$$
\frac{0.141 - 0.109}{0.141} = 22.69\%
$$
Total loss is
$$
\frac{1.187 - 2 * 0.095}{1.187} = 92\%
$$
This is called by the bending loss caused by the bendings in the waveguid. Here owing to the 90 degree bending, the losses are huge which are obuiously not we want.

Here we can see how the wave propagate and the leakage.


### Y-Splitter
To reduce the bending loss, a new scheme is than ideated.
It consist of two oblique waveguides and is connected to the main waveguide by a waveguide taper.

The resulting $power/m^2$ is

1.319 0.444 0.448

The splitt ratios are$\frac{0.444}{0.444 + 0.448} = 49.8\%$ and $50.2\%$. The total loss is
$$
\frac{1.319 - 0.444 - 0.448}{1.319} = 32.4\%
$$
The loss is still high due to the mismatch of waveguide taper and the bending loss.

### YS-Splitter
[Syahriar mentioned that](10.1109/RFM.2008.4897467) S-shaped wave guide can reduce the bending loss due to its continuous first-derivatice feature. Here I try to examine this idea.
The resulting $power/m^2$ is

1.411 0.685 0.698

The splitt ratios are$\frac{0.685}{0.685 + 0.698} = 49.5\%$ and $50.5\%$.
The total loss is
$$
\frac{1.411 - 0.685 - 0.698}{1.411} = 1.98\%
$$

The loss has significant improvement.
(Note that this is a worser resolution version due to the upload constraint of gif.
The the dispersion is due to geometry defects which can be sorted by increasing resolution)

### Coupling-Splitter
Beside traditional-designed splitter, this simple demonstration of splitter using couple-mode theory is implement and it shows this thought is feasible.

This scheme consist of two adjacent straight wave guides used to transfer the guided field from one waveguide to another. The transfered field is then splited by a S-shape waveguide. By changing the length of adjacent part, we can adjust the split ratio.

The resulting $power/m^2$ is

1.404 0.767 0.63

The splitt ratios are$\frac{0.767}{0.767 + 0.63} = 55\%$ and $45\%$.

The total loss is
$$
\frac{1.404 - 0.767 - 0.63}{1.404} = 0.5\%
$$
The loss is small, and by properly adjust the length of adjacent region, it may reaches 50-50 splitt ratios.
### MMI Splitter
Another idea is to split the light with Multi mode interferometer (MMI). By using the [talbot effect](https://www.sciencedirect.com/science/article/pii/B9780125250962500039), a $m \times m$ coupler can than be designed.

The MMI device is made of two sections of single mode waveguide, and one section of multi mode wave guide with width $W_{multimode} = 10 W_{singlemode}$.
The resulting $power/m^2$ is

1.333 0.675 0.624

The splitt ratios are$\frac{0.675}{0.675 + 0.624} = 52\%$ and $48\%$.

The total loss is
$$
\frac{1.333 - 0.675 - 0.624}{1.404} = 3.4\%
$$

## Conclusion

In this report we go through some fundemental knowledge of waveguide and FDTD. In the later section we investigate and simulate through different types of 50-50 splitter. Giving a brief inspect of PIC components.