[[ coldCLOUD ]]

# 2012.10 3D Waveguide Model

In integrated optics it's often very crucial to know the electromagnetic (EM) field distribution inside and around the optical waveguide one is considering.

This may be done by defining the waveguide geometry in terms of the refractive index $n$ and its spatial and spectral dependence.

Now, the temporal evolution of a given EM field source may be found be solving Maxwell's equations for that given structure. There are several different methods for solving Maxwell's equations numerically and these are used by different numerical solvers.

One of these methods is the Finite-Difference Time-Domain (FDTD) method and it's used by one (freely available) numerical solver called Meep.

## Model

The simple model used here is a rectangular waveguide of Silicon Nitride (${\rm SiN}, \; n = 1.98$) buried in Silicon Dioxide (${\rm SiO_2}, \; n= 1.46$).

In terms of the dimensionless length-scale parameter $a$ in meep the waveguide crosssection is $w = 0.5, \; h = 0.25, \; l = 5$.

The structure is pumped using a monochromatic continuous-wave dipole source of frequency $f = 0.850$ that only has an electric field component $E_y$ along the y-axis - i.e. a transverse electric (TE) field.

The implementation in Meep is done using Scheme and looks like this,

; Materials
(define SiO2 (make dielectric (index 1.46)))
(define SiN (make dielectric (index 1.98)))

; Waveguide
(define-param mwg SiN) ; waveguide material
(define-param wwg 0.500) ; waveguide width
(define-param hwg 0.250) ; waveguide height
(define-param lwg 5.000) ; initial waveguide length

(define-param mcl SiO2) ; cladding material

; Light
(define-param lfc 0.850) ; light pulse center frequency
(define-param ldf 0.0)  ; light pulse width [in frequency]

(define-param tpml 1) ; thickness of PML

; Resolution
(define-param res 20) ; resolution of computational cell

; Default material
(set! default-material mcl)

; Run with sub-pixel averaging?
(set! eps-averaging? false)

; Calculation of parameters
(define xsz (+ lwg (* 2 xpad) (* 2 tpml)) ) ; x-size of computational cell
(define ysz (+ wwg (* 2 ypad) (* 2 tpml)) ) ; y-size of computational cell
(define zsz (+ hwg (* 2 zpad) (* 2 tpml)) ) ; z-size of computational cell

; SET COMPUTATIONAL CELL SIZE
(set! geometry-lattice
(make lattice (size xsz ysz zsz))
)

; MAKE GEOMETRY
(set! geometry
(list
; Build WG and extend it thru PMLs and PADs
(make block (center 0 0 0) (size xsz wwg hwg) (material mwg))
)
)

; SET PMLs
(set! pml-layers (list (make pml (thickness tpml))))

; SET RESOLUTION
(set-param! resolution res)

; SET SOURCES
(set! sources
(list
(make source
(src (make continuous-src (frequency lfc) (end-time 200) ))
(component Ey) (size 0 0 0) (center (- (+ tpml (* 0.1 xpad)) (/ xsz 2) ) 0 0)
)
)
)

(run-until 60
(at-beginning output-epsilon)
(at-every 10 output-efield-y)
)

A block-by-block explanation of the model follows here,

• definition of material refractive indices
• waveguide geometry
• cladding material (here used to set the default material later)
• properties of light source
• define padding zone and perfectly matched layers (PML) (to prevent reflection from computational cell (CC) boundaries)
• CC resolution (number of points per length-scale $a$)
• disable sub-pixel averaging (sharp index differences are smoothed out - not desirable here)
• calculate the size of CC given the geometry, padding zones and PML's
• set CC size
• set geometry
• set PML's
• set CC resolution
• set light source (position, type, component, etc)
• set run time and what to export (here refractive index structure and the relevant EM field component)

Note: if you run this code be careful not to make the CC too big or the resolution too high - it will slow down your system. It's already computational heavy and may be speeded up by running the code in parallel.

## Visualization

The exported files are stored in a standard data file format called HDF5 and it's supported by many other software platforms, e.g. Mathematica.

Using the electric field component data file at time step $t = 40$ we have the following situation.

The electric field component intensity in different cross sections is shown below,

Green lines mark the position of the waveguide edges, while the red and black lines on the contour plot mark the longitudinal and transverse field intensity distributions shown in the two other plots. The axes on the plots are pixel - or data point - number.

Visualizing this information in a 3D fashion we have,

### Mathematica Code

Here's a small note on the Mathematica code used to combine the 2D contours plots into a 3D object.

Using these three functions we may convert a 2D plot into a cross section of a 3D plot,

(* Insert a 2D plot into the xy, xz, or yz plane of a 3D object - x0, y0, and z0 are offsets, while opacity is the opacity of the 2D plot*)
Make3Dxy[plot_, {x0_, y0_, z0_}, opacity_] := Module[{newplot},
newplot = First@Graphics[plot];
newplot = N@newplot /. {x_?AtomQ, y_?AtomQ} -> {x + x0, y + y0, z0};
newplot /.
GraphicsComplex[xx__] -> {Opacity[opacity], GraphicsComplex[xx]}
];
Make3Dxz[plot_, {x0_, y0_, z0_}, opacity_] := Module[{newplot},
newplot = First@Graphics[plot];
newplot = N@newplot /. {x_?AtomQ, z_?AtomQ} -> {x + x0, y0, z + z0};
newplot /.
GraphicsComplex[xx__] -> {Opacity[opacity], GraphicsComplex[xx]}
];
Make3Dyz[plot_, {x0_, y0_, z0_}, opacity_] := Module[{newplot},
newplot = First@Graphics[plot];
newplot = N@newplot /. {y_?AtomQ, z_?AtomQ} -> {x0, y + y0, z + z0};
newplot /.
GraphicsComplex[xx__] -> {Opacity[opacity], GraphicsComplex[xx]}
];

This code is based on the idea described here.

Here I've just extended it to a insert 2D plot into any Cartesian plane.

What the code does is to,

• open up the 2D graphics object
• look for the list of coordinate points, maniplulate them and replace the list by a 3D coordinate list
• recombine it into a graphics object while setting the opacity 