Sitemap Recent changes

This is an old revision of the document!

Plane wave expansion

Eigenmode decomposition is very popular tool to analyze electromagnetic ﬁelds inside arbitrary optical systems. In the absence of nonlinearity, any electromagnetic propagation could be represented as a superposition of eigenmodes.

One approach of eigenmodes calculation relies on expanding the eigenmode in series of specially chosen basic functions. Such decomposing leads to finite matrix eigenproblem which can be solved using standard linear algebra.

Different choices of basic functions can be applied correspondingly to considered problem. For example, Chebyshev or Legendre polynomials can be used for finite optical systems, while Fourier series can be used for infinite photonic crystals. Second choice is a core of plane wave expansion (PWE) method that is very suitable for calculation of photonic crystals band structures.

We implemented plane wave expansion (PWE) method as a part of Electromagnetic Template Library (EMTL). Note that EMTL was initially developed mostly for FDTD simulations. Contrary to FDTD, our PWE implementation is much more restricted: currently it does not support calculation of anisotropic or dispersive materials, it is not parallel, etc. However, it can be used for efficient simulation of 2D and 3D photonic crystals.

Below we present illustrative tutorial how to use EMTL for PWE simulations. There is some common syntax with that one used for FDTD, especially concerning how to construct geometry objects, media and media objects.

Bandstructure of 2D photonic crystals

Here we show how to calculate bandstructure of 2D photonic crystal.

To use EMTL for PWE simulations you should include file plane_wave_exp.h to your code

#include "plane_wave_exp.h"

EMTL should be initialized with command line arguments

int main(int argc, char **argv){
emInit(argc,argv);
return 0;
}
Function emInit has 3rd optional arguments. You can specify directory for output text files there. This argument is empty strings by default (this means that all files are recorded to directory with your executable file).

The main class for numerical experiment is pwExperiment declared in photonic/plane_wave_exp.h. Lets make object of this class.

pwExperiment task;

We will calculate band structure for the TM waves inside a square lattice of circular rods (epsilon=11.9, radius 0.2) in air background.

We define periodic lattice and add circular rod at the central position of the periodic cell:

task.SetLattice(Vector_3(1,0,0),Vector_3(0,1,0));
task.AddObject(emMedium(11.9), GetCylinder(Vector_3(.5,.5,0),Vector_3(0,0,1),.2));

We specify TM polarization (it makes sense only for 2D):

task.SetPolarization("TM");

We are interested in first 10 bands:

task.SetBandsNumber(10);

We will use 50 spatial mesh steps along each dimension. Spatial mesh is used to calculate Fourier representation of inverse dielectric permittivity.

task.SetMeshStepsNumber(50);

We will use plane waves along each dimension (5 for G>0, 5 for G<0, and one is G=0). Total number of plane waves is .

task.SetPlaneWavesNumber(5);

This is our square lattice and first Brillouin zone.

We calculate eigenfrequencies for a number of wavevectors on the path .

int rsl=10;
for(int i=0;i<rsl;i++)
for(int i=0;i<rsl;i++)
for(int i=rsl;i>=0;i--)
task.Calculate();

As a result of calculation we get tabular file 'eig.d' with first 10 bands. Table in this file has following columns: number of wavector, three coordinates of wavevector, band frequencies. Open this file by any text editor and see how is it organized.

You can plot bandstructure from this file using popular gnuplot program that can be downloaded from http://www.gnuplot.info.

For example, to plot first 10 bands you should open gnupolt, move to directory with eig.d file using cd command (specify directory name in quotes) and type

gnuplot> unset key
gnuplot> set ylabel 'a/lambda'
gnuplot> set xtics ("G" 0, "X" 10, "M" 20, "G" 30)
gnuplot> p 'eig.d' u 0:5 w l, 'eig.d' u 0:6 w l, 'eig.d' u 0:7 w l, 'eig.d' u 0:8 w l, 'eig.d' u 0:9 w l, 'eig.d' u 0:10 w l, 'eig.d' u 0:11 w l, 'eig.d' u 0:12 w l, 'eig.d' u 0:13 w l, 'eig.d' u 0:14 w l

This is what you will get:

We hide legend of the plot using 'unset key' gnuplot comand. We also label wavevector numbers 0, 10, 20, 30 by corresponding letters.

In the file bandgaps.d you will find bandgap frequency ranges (if they exist) corresponding to number of specified bands.

Waveguide modes. Supercell technique.

Translation symmetry of photonic crystal can be broken by a defect. Defects may permit localized modes to exist, with frequencies inside photonic band gaps. Using various types of defects one can localize light in the vicinity of a single point (an optical cavity), direct along a particular way (a waveguide), or trap it at a 2D surface. In this section we will illustrate this using example of waveguide formed by linear chain of defect rods inside square lattice of rods from previous section.

To calculate band structures for a photonic crystal with a linear defect we use supercell technique. We assume a periodic structure with a large primitive cell (called a supercell), which includes more than one unit cell. The supercell contains a defect rod (radius chosen differently from the other rods). This supercell is then repeated periodically in a rectangular lattice to form an array of linear waveguides (rows with defect rods) separated by many unit cells of photonic crystal. The effect of coupling between two adjacent waveguides can be made negligible by using a very large supercell.

This is example of applying supercell technique

#include "plane_wave_exp.h"

int main(int argc,char **argv){

emInit(argc,argv);

int super=4;
for(int i=0;i<super;i++){
}

int rsl=10;
for(int i=0;i<=rsl;i++)

return 0;
}

Supercell consists on 2*super+1 unit cells. Defect cell is placed at the center of the supercell. Using this code you can calculate the eigenfrequencies for a number of wave vectors on the path to .

This is the result for first 15 bands (supercell=4).

gnuplot> unset key
gnuplot> set xlabel 'ky, 2pi/a'
gnuplot> set ylabel 'a/lambda'
gnuplot> set xtics
gnuplot> p 'eig.d' u 3:5 w l, 'eig.d' u 3:6 w l, 'eig.d' u 3:7 w l, 'eig.d' u 3:8 w l, 'eig.d' u 3:9 w l, 'eig.d' u 3:10 w l, 'eig.d' u 3:11 w l, 'eig.d' u 3:12 w l, 'eig.d' u 3:13 w l, 'eig.d' u 3:14 w l, 'eig.d' u 3:15 w l, 'eig.d' u 3:16 w l, 'eig.d' u 3:17 w l, 'eig.d' u 3:18 w l, 'eig.d' u 3:19 w l