- MEEP Project #1 — Light-Extraction Efficiency of Organic Light Emitting Diodes (OLEDs)
- MPB Project #1 — Modes of Silicon on Insulator (SOI) Strip Waveguides
- MEEP Project #2 — Optimizing Far-Field Intensity of SOI Bragg Grating Outcouplers
- MPB Project #2 — Band Gap of Photonic-Crystal Nanobeam Waveguide
- MEEP Project #3 — Resonant Modes of a Photonic-Crystal Nanobeam Cavity

There are three key features involved in developing an accurate model. (1)

` ````
(set-param! resolution 100) ; pixels/unit length (1 μm)
```

We specify the frequency bounds of the Gaussian pulse using its minimum and maximum wavelengths in vacuum.
` ````
(define-param lambda-min 0.4) ; minimum source wavelength
(define-param lambda-max 0.8) ; maximum source wavelength
(define fmin (/ lambda-max)) ; minimum source frequency
(define fmax (/ lambda-min)) ; maximum source frequency
(define fcen (* 0.5 (+ fmin fmax))) ; source frequency center
(define df (- fmax fmin)) ; source frequency width
```

The thickness of each layer is a parameter. This is useful if we want to optimize the design. We use typical values for OLEDs. The length of the absorbing layer should be at least half the largest wavelength in the simulation to ensure negligible reflections. This is because incident waves are absorbed twice through a round-trip reflection from the hard-wall boundaries of the computational cell.
` ````
(define-param tABS lambda-max) ; absorber/PML thickness
(define-param tGLS 1) ; glass thickness
(define-param tITO 0.1) ; indium tin oxide thickness
(define-param tORG 0.1) ; organic thickness
(define-param tAl 0.1) ; aluminum thickness
; length of computational cell along Z
(define sz (+ tABS tGLS tITO tORG tAl))
```

Since this is a 3d simulation, we also need to specify the length of the computational cell in the transverse directions X and Y. This is the length of the high-index organic/ITO waveguide. We choose a value of several wavelengths. The choice of the waveguide length has a direct impact on the results.
` ````
; length of non-absorbing region of computational cell in X and Y
(define-param L 4)
(set! geometry-lattice (make lattice (size (+ L (* 2 tABS)) (+ L (* 2 tABS)) sz)))
```

Overlapping the lossy-metal aluminum with a perfectly-matched layer (PML) sometimes leads to field instabilities due to backward-wave SPPs. MEEP provides an alternative absorber which tends to be more stable. We use an absorber in the X and Y directions and a PML for the outgoing waves in the glass substrate. The metal cathode on top of the device either reflects or absorbs the incident light. No light is transmitted. Thus, we need a PML along just one side in the Z direction.
` ````
(set! pml-layers (list
(make absorber (thickness tABS) (direction X))
(make absorber (thickness tABS) (direction Y))
(make pml (thickness tABS) (direction Z) (side High))))
```

Next, we define the material properties and set up the geometry consisting of the four-layer stack.
` ````
(define-param nORG 1.75)
(define ORG (make medium (index nORG)))
(define-param nITO 1.8)
(define ITO (make medium (index nITO)))
(define-param nGLS 1.45)
(define GLS (make medium (index nGLS)))
; conversion factor for eV to 1/um
(define eV-um-scale (/ 1.23984193))
; Al, from Rakic et al., Applied Optics, vol. 32, p. 5274 (1998)
(define Al-eps-inf 1)
(define Al-plasma-frq (* 14.98 eV-um-scale))
(define Al-f0 0.523)
(define Al-frq0 1e-10)
(define Al-gam0 (* 0.047 eV-um-scale))
(define Al-sig0 (/ (* Al-f0 (sqr Al-plasma-frq)) (sqr Al-frq0)))
(define Al-f1 0.050)
(define Al-frq1 (* 1.544 eV-um-scale)) ; 803 nm
(define Al-gam1 (* 0.312 eV-um-scale))
(define Al-sig1 (/ (* Al-f1 (sqr Al-plasma-frq)) (sqr Al-frq1)))
(define Al-f2 0.166)
(define Al-frq2 (* 1.808 eV-um-scale)) ; 686 nm
(define Al-gam2 (* 1.351 eV-um-scale))
(define Al-sig2 (/ (* Al-f2 (sqr Al-plasma-frq)) (sqr Al-frq2)))
(define Al-f3 0.030)
(define Al-frq3 (* 3.473 eV-um-scale)) ; 357 nm
(define Al-gam3 (* 3.382 eV-um-scale))
(define Al-sig3 (/ (* Al-f3 (sqr Al-plasma-frq)) (sqr Al-frq3)))
(define Al (make medium (epsilon Al-eps-inf)
(E-polarizations
(make drude-susceptibility
(frequency Al-frq0) (gamma Al-gam0) (sigma Al-sig0))
(make lorentzian-susceptibility
(frequency Al-frq1) (gamma Al-gam1) (sigma Al-sig1))
(make lorentzian-susceptibility
(frequency Al-frq2) (gamma Al-gam2) (sigma Al-sig2))
(make lorentzian-susceptibility
(frequency Al-frq3) (gamma Al-gam3) (sigma Al-sig3)))))
(set! geometry (list
(make block (material GLS) (size infinity infinity (+ tABS tGLS))
(center 0 0 (- (* 0.5 sz) (* 0.5 (+ tABS tGLS)))))
(make block (material ITO) (size infinity infinity tITO)
(center 0 0 (- (* 0.5 sz) tABS tGLS (* 0.5 tITO))))
(make block (material ORG) (size infinity infinity tORG)
(center 0 0 (- (* 0.5 sz) tABS tGLS tITO (* 0.5 tORG))))
(make block (material Al) (size infinity infinity tAl)
(center 0 0 (- (* 0.5 sz) tABS tGLS tITO tORG (* 0.5 tAl))))))
```

We use a set of point-dipole sources, each with random phase but fixed polarization, distributed along a line within the middle of the organic layer. The number of point sources is a parameter with a default of 10. The polarization of the source has an important effect on the results which we will investigate. Instead of using a random polarization, which would require multiple runs to obtain a statistical average, the sources can be separated into components which are parallel (directions X and Y) and perpendicular (Z) to the layers. Thus, we will need to do just two sets of simulations. These are averaged to obtain results for random polarization.
` ````
; random number generator: uniformly distributed in [-1,1]
(define random-num (lambda ()
(let ((time (gettimeofday)))
(set! *random-state* (seed->random-state (+ (car time) (cdr time)))))
(if (> (random:uniform) 0.5) (random:uniform) (* -1 (random:uniform)))))
(define-param src-cmpt Ex) ; current source component
(define-param num-src 10) ; number of point sources
(set! sources
(map (lambda (cz)
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component src-cmpt) (center 0 0 (- (* 0.5 sz) tABS tGLS tITO (* 0.4 tORG) (* 0.2 cz tORG)))
(amplitude (exp (* 0+2i pi (abs (random-num)))))))
(arith-sequence (/ num-src) (/ num-src) num-src)))
```

We can also exploit the two mirror symmetry planes of the sources and the structure in order to reduce the size of the computation by a factor four. The phase of the mirror symmetry planes depends on the polarization of the source. Three separate cases are necessary.
` ````
(if (= src-cmpt Ex)
(set! symmetries (list
(make mirror-sym (direction X) (phase -1))
(make mirror-sym (direction Y) (phase +1)))))
(if (= src-cmpt Ey)
(set! symmetries (list
(make mirror-sym (direction X) (phase +1))
(make mirror-sym (direction Y) (phase -1)))))
(if (= src-cmpt Ez)
(set! symmetries (list
(make mirror-sym (direction X) (phase +1))
(make mirror-sym (direction Y) (phase +1)))))
```

Three sets of flux monitors are required to capture the power in the sources, glass substrate, and high-index waveguide. We place a set of six flux planes along the sides of a box to enclose the sources. This is used to compute the total power in the device. These monitors must be placed entirely within the organic layer. One flux plane is required to compute the power in the glass substrate. Four flux planes are required to capture the total power of the modes in the high-index waveguide forming the organic/ITO layers. These modes are delocalized and extend slightly beyond these two layers. This requires the height of the flux planes to be larger than the combined thickness of the layers as shown in the figure above. Also, any component of the guided mode which extends into the aluminum will be absorbed due to screening effects of the charges in the metal. The longer the waveguide is made, the more these guided modes will be absorbed. This is why the length of the waveguide (L) should not be made so large that the waveguide component of the total power is zero.
` ````
; number of frequency bins for DFT fields
(define-param nfreq 50)
; surround source with a six-sided box of flux planes
(define srcbox-width 0.05)
(define srcbox-top (add-flux fcen df nfreq
(make flux-region (size srcbox-width srcbox-width 0) (direction Z)
(center 0 0 (- (* 0.5 sz) tABS tGLS)) (weight +1))))
(define srcbox-bot (add-flux fcen df nfreq
(make flux-region (size srcbox-width srcbox-width 0) (direction Z)
(center 0 0 (- (* 0.5 sz) tABS tGLS tITO (* 0.8 tORG))) (weight -1))))
(define srcbox-xp (add-flux fcen df nfreq
(make flux-region (size 0 srcbox-width (+ tITO (* 0.8 tORG))) (direction X)
(center (* 0.5 srcbox-width) 0 (- (* 0.5 sz) tABS tGLS (* 0.5 (+ tITO (* 0.8 tORG))))) (weight +1))))
(define srcbox-xm (add-flux fcen df nfreq
(make flux-region (size 0 srcbox-width (+ tITO (* 0.8 tORG))) (direction X)
(center (* -0.5 srcbox-width) 0 (- (* 0.5 sz) tABS tGLS (* 0.5 (+ tITO (* 0.8 tORG))))) (weight -1))))
(define srcbox-yp (add-flux fcen df nfreq
(make flux-region (size srcbox-width 0 (+ tITO (* 0.8 tORG))) (direction Y)
(center 0 (* 0.5 srcbox-width) (- (* 0.5 sz) tABS tGLS (* 0.5 (+ tITO (* 0.8 tORG))))) (weight +1))))
(define srcbox-ym (add-flux fcen df nfreq
(make flux-region (size srcbox-width 0 (+ tITO (* 0.8 tORG))) (direction Y)
(center 0 (* -0.5 srcbox-width) (- (* 0.5 sz) tABS tGLS (* 0.5 (+ tITO (* 0.8 tORG))))) (weight -1))))
; padding for flux box to capture extended waveguide mode
(define fluxbox-dpad 0.05)
; upward flux into glass substrate
(define glass-flux (add-flux fcen df nfreq
(make flux-region (size L L 0) (direction Z)
(center 0 0 (- (* 0.5 sz) tABS (- tGLS fluxbox-dpad))) (weight +1))))
; surround ORG/ITO waveguide with four-sided box of flux planes
; NOTE: waveguide mode extends partially into Al cathode and glass substrate
(define wvgbox-xp (add-flux fcen df nfreq
(make flux-region (size 0 L (+ fluxbox-dpad tITO tORG fluxbox-dpad)) (direction X)
(center (* 0.5 L) 0 (- (* 0.5 sz) tABS tGLS (* 0.5 (+ tITO tORG)))) (weight +1))))
(define wvgbox-xm (add-flux fcen df nfreq
(make flux-region (size 0 L (+ fluxbox-dpad tITO tORG fluxbox-dpad)) (direction X)
(center (* -0.5 L) 0 (- (* 0.5 sz) tABS tGLS (* 0.5 (+ tITO tORG)))) (weight -1))))
(define wvgbox-yp (add-flux fcen df nfreq
(make flux-region (size L 0 (+ fluxbox-dpad tITO tORG fluxbox-dpad)) (direction Y)
(center 0 (* 0.5 L) (- (* 0.5 sz) tABS tGLS (* 0.5 (+ tITO tORG)))) (weight +1))))
(define wvgbox-ym (add-flux fcen df nfreq
(make flux-region (size L 0 (+ fluxbox-dpad tITO tORG fluxbox-dpad)) (direction Y)
(center 0 (* -0.5 L) (- (* 0.5 sz) tABS tGLS (* 0.5 (+ tITO tORG)))) (weight -1))))
```

Finally, with the structure and sources set up, we run the simulation until the fields in the device have sufficiently decayed away. Afterwards, the flux data is written to standard output as a series of columns.
` ````
(run-sources+ (stop-when-fields-decayed 50 src-cmpt (vector3 0 0 (- (* 0.5 sz) tABS tGLS tITO (* 0.5 tORG))) 1e-8))
(display-fluxes srcbox-top srcbox-bot srcbox-xp srcbox-xm srcbox-yp srcbox-ym glass-flux wvgbox-xp wvgbox-xm wvgbox-yp wvgbox-ym)
```

We use the following Bash script to launch a parallel MEEP simulation using 8 processors.
` ````
#!/bin/bash
src_cmpt=$1;
L=$2;
mpirun -np 8 meep-mpi src-cmpt=E${src_cmpt} L=${L} oled-ext-eff.ctl > oled-flux.out;
grep flux1: oled-flux.out |cut -d , -f2- > oled-flux.dat;
```

This script is executed from the shell terminal using a Jx source and waveguide length L of 4.
` ````
chmod u+x run_oled.sh
nohup ./run_OLED.sh x 4 &> /dev/null &
```

This takes about 3 hours to run on an Intel Xeon 2.60 GHz machine. Once the simulation is complete, we plot its output using a Jupyter/IPython notebook. The total power in each device component is the sum of the relevant columns of the output.
` ````
from numpy import *
from scipy import interpolate
import matplotlib.pyplot as plt
f = genfromtxt("oled-flux.dat", delimiter=",")
total_power = sum(f[:,1:6],axis=1)
glass_power = f[:,7]
waveguide_power = sum(f[:,8:11],axis=1)
glass = glass_power/total_power
waveguide = waveguide_power/total_power
aluminum = 1-glass-waveguide
lambdas_linear = linspace(0.4,0.8,100)
lambdas = 1/f[:,0]
g_linear = interpolate.interp1d(lambdas,glass,kind='cubic')
w_linear = interpolate.interp1d(lambdas,waveguide,kind='cubic')
a_linear = interpolate.interp1d(lambdas,aluminum,kind='cubic')
glass_linear = g_linear(lambdas_linear)
waveguide_linear = w_linear(lambdas_linear)
aluminum_linear = a_linear(lambdas_linear)
plt.plot(lambdas_linear,glass_linear,'b-',label='glass');
plt.plot(lambdas_linear,waveguide_linear,'r-',label='organic + ITO');
plt.plot(lambdas_linear,aluminum_linear,'g-',label='aluminum');
plt.xlabel("wavelength (um)"); plt.ylabel("fraction of total power");
plt.axis([0.4, 0.8, 0, 1]);
plt.xticks([t for t in arange(0.4,0.9,0.1)]);
plt.legend(loc='upper right');
plt.show()
print("glass: %0.6f (mean), %0.6f (std. dev.)" % (mean(glass_linear),std(glass_linear)))
print("waveguide: %0.6f (mean), %0.6f (std. dev.)" % (mean(waveguide_linear),std(waveguide_linear)))
print("aluminum: %0.6f (mean), %0.6f (std. dev.)" % (mean(aluminum_linear),std(aluminum_linear)))
```

The left figure shows the device structure. The silicon waveguide has a rectangular cross section with width

` ````
(set-param! resolution 64) ; pixels/unit length (1 μm)
```

The width and height of the waveguide are design parameters. Most SOI wafers have a 220 nm thick silicon layer which is used as the default value for the height. The default width is 500 nm.
` ````
(define-param w 0.50) ; waveguide width
(define-param h 0.22) ; waveguide height
```

This is a 3d calculation based on a 2d computational cell in Y and Z. The size of the computational cell needs to be large enough such that the fields are negligible at the edge of the computational cell. This ensures that MPB's periodic boundaries do not affect the results. Since the guided mode is localized within the silicon, the fields outside are exponentially decaying. There is a tradeoff to increasing the size of the computational cell: the size of the calculation also increases.
` ````
(define-param sc-y 2) ; supercell width
(define-param sc-z 2) ; supercell height
(set! geometry-lattice (make lattice (size no-size sc-y sc-z)))
```

Next, we specify the material properties of silicon and silicon dioxide. Silicon is lossless at 1.55 μm which is an important feature for these kinds of applications. The large refractive index contrast with silicon dioxide also produces strong mode confinement. Typical values for the refractive indicies are used as default.
` ````
(define-param nSi 3.45)
(define-param nSiO2 1.45)
(define Si (make dielectric (index nSi)))
(define SiO2 (make dielectric (index nSiO2)))
```

Since the oxide layer in which the waveguide modes are exponentially decaying is several microns in thickness, the fields are negligible at the silicon substrate. Thus, the substrate can be omitted from the device geometry. The geometry consists of just two objects: (1) a rectangular silicon waveguide at the center of the computational cell and (2) a planar oxide layer which fills the bottom half.
` ````
(set! geometry (list
(make block (size infinity w h)
(center 0 0 0) (material Si))
(make block (size infinity infinity (* 0.5 (- sc-z h)))
(center 0 0 (* 0.25 (+ sc-z h))) (material SiO2))))
```

The first four bands will be computed even though the single-mode property can be determined from only the first and second bands. The range of propagation wavevectors over which to compute the bands needs to be sufficiently large to contain modes at 1.55 μm. Finding suitable values typically requires a bit of trial and error.
` ````
(set-param! num-bands 4)
(define-param num-k 20)
(define-param k-min 0.1)
(define-param k-max 2.0)
(set! k-points (interpolate num-k (list (vector3 k-min) (vector3 k-max))))
```

We can now call the `run`

function to compute the bands. The parity information in the Y direction for each mode is also output. This will help us identify properties of the modes since the fundamental mode has even symmetry.
` ````
(run display-yparities)
```

After the bands have been computed, we identify the mode(s) at 1.55 μm. This involves an inverse calculation where given the mode frequency (in vacuum), we compute its wavevector using MPB's built-in `find-k`

routine. The Poynting vector along X of the fields is also output as an HDF5 file.
` ````
(define-param f-mode (/ 1.55)) ; frequency corresponding to 1.55 um
(define-param band-min 1)
(define-param band-max 1)
(define-param kdir (vector3 1 0 0))
(define-param tol 1e-6)
(define-param kmag-guess (* f-mode nSi))
(define-param kmag-min (* f-mode 0.1))
(define-param kmag-max (* f-mode 4))
(find-k ODD-Y f-mode band-min band-max kdir tol
kmag-guess kmag-min kmag-max
output-poynting-x display-group-velocities)
```

This shell script is used to launch the simulation which takes a few seconds using a single Intel Xeon 2.60 GHz core. Afterwards, `h5topng`

, part of the h5utils package, is used to generate a PNG image from the raw simulation data stored in the HDF5 file.
` ````
#!/bin/bash
mpb strip-wvg.ctl > strip-wvg-bands.out;
grep freqs: strip-wvg-bands.out |cut -d , -f3,7- |sed 1d > strip-wvg-bands.dat;
h5topng -o wvg_power.png -x 0 -d x.r -vZc bluered -C strip-wvg-epsilon.h5 strip-wvg-flux.v.k01.b01.x.yodd.h5;
```

Finally, we plot the results from the output and add the light line for the oxide.
` ````
from numpy import *
import matplotlib.pyplot as plt
f = genfromtxt("strip-wvg-bands.dat", delimiter=",")
plt.plot(f[:,0],f[:,1:],'b-',f[:,0],f[:,0]/1.45,'k-');
plt.xlabel("wavevector k_x (units of 2\pi\mum^{-1})");
plt.ylabel("frequency (units of 2\pix30^{14} Hz)");
plt.axis([0,2,0,1]);
plt.show()
```

This plot is shown in the right figure above. There is a single mode at 1.55 μm for the lowest band indicated with the dotted green line. The cutoff frequency of the second band lies above this mode frequency. The inset shows the Poynting vector along X (S`find-k`

routine.The outcoupler design is based on a concentric Bragg grating with angled sides, shown in the figures below. The input port is an SOI strip waveguide which is connected to the Bragg grating.

The figure below shows the device cross section in the XY plane of the computational cell. There are two parameters used to design the Bragg grating outcoupler: periodicity

The simulation script consists of the device geometry, source, near-field surface (used to compute the far fields), time-stepping routine, and far-field output. First, the device geometry which includes the Bragg grating, waveguide, and substrate is defined using a combination of various shape objects. The entire structure is parameterized though only two parameters are optimized in this example.

` ````
(set-param! resolution 20) ; pixels/unit length (1 um)
(define-param h 0.22) ; waveguide height
(define-param w 0.5) ; waveguide width
(define-param a 1.0) ; Bragg grating periodicity/lattice parameter
(define-param d 0.5) ; Bragg grating thickness
(define-param N 5) ; number of grating periods
(set! N (+ N 1))
(define nSi 3.45)
(define Si (make medium (index nSi)))
(define nSiO2 1.45)
(define SiO2 (make medium (index nSiO2)))
(define-param sxy 16)
(define-param sz 4)
(set! geometry-lattice (make lattice (size sxy sxy sz)))
; rings of Bragg grating
(set! geometry (append geometry
(map (lambda (n)
(list
(make cylinder (material Si) (center 0 0 0) (radius (* n a)) (height h))
(make cylinder (material air) (center 0 0 0) (radius (- (* n a) d)) (height h))))
(arith-sequence N -1 N))))
(set! geometry (apply append geometry))
; remove left half of Bragg grating rings to form semi circle
(set! geometry (append geometry (list
(make block (material air) (center (* -0.5 (* N a)) 0 0) (size (* N a) (* 2 N a) h))
(make cylinder (material Si) (center 0 0 0) (radius (- a d)) (height h)))))
; angle sides of Bragg grating
; rotation angle of sides relative to Y axis (degrees)
(define-param rot-theta 0)
(set! rot-theta (deg->rad (- rot-theta)))
(define pvec (vector3 0 (* 0.5 w) 0))
(define cvec (vector3 (* -0.5 N a) (+ (* 0.5 N a) (* 0.5 w)) 0))
(define rvec (vector3- cvec pvec))
(define rrvec (rotate-vector3 (vector3 0 0 1) rot-theta rvec))
(set! geometry (append geometry (list (make block
(material air)
(center (vector3+ pvec rrvec)) (size (* N a) (* N a) h)
(e1 (rotate-vector3 (vector3 0 0 1) rot-theta (vector3 1 0 0)))
(e2 (rotate-vector3 (vector3 0 0 1) rot-theta (vector3 0 1 0)))
(e3 (vector3 0 0 1))))))
(set! pvec (vector3 0 (* -0.5 w) 0))
(set! cvec (vector3 (* -0.5 N a) (- (+ (* 0.5 N a) (* 0.5 w))) 0))
(set! rvec (vector3- cvec pvec))
(set! rrvec (rotate-vector3 (vector3 0 0 1) (- rot-theta) rvec))
(set! geometry (append geometry (list (make block
(material air)
(center (vector3+ pvec rrvec)) (size (* N a) (* N a) h)
(e1 (rotate-vector3 (vector3 0 0 1) (- rot-theta) (vector3 1 0 0)))
(e2 (rotate-vector3 (vector3 0 0 1) (- rot-theta) (vector3 0 1 0)))
(e3 (vector3 0 0 1))))))
; input waveguide
(set! geometry (append geometry (list
(make block (material air) (center (* -0.25 sxy) (+ (* 0.5 w) (* 0.5 a)) 0) (size (* 0.5 sxy) a h))
(make block (material air) (center (* -0.25 sxy) (- (+ (* 0.5 w) (* 0.5 a))) 0) (size (* 0.5 sxy) a h))
(make block (material Si) (center (* -0.25 sxy) 0 0) (size (* 0.5 sxy) w h)))))
; substrate
(set! geometry (append geometry (list (make block
(material SiO2)
(center 0 0 (+ (* -0.5 sz) (* 0.25 (- sz h))))
(size infinity infinity (* 0.5 (- sz h)))))))
; surround the entire computational cell with PML
(define-param dpml 1.0)
(set! pml-layers (list (make pml (thickness dpml))))
```

Next, the eigenmode source is defined which is based on calling MPB to compute a definite-frequency mode and using it as the amplitude profile. We use a Gaussian pulsed source, with center wavelength corresponding to 1.55 μm, such that the fields eventually decay away due to absorption by the PMLs and the simulation can be terminated.
` ````
; mode frequency
(define-param fcen (/ 1.55))
(set! sources (list (make eigenmode-source
(src (make gaussian-src (frequency fcen) (fwidth (* 0.2 fcen))))
(component Ey)
(size 0 (- sxy (* 2 dpml)) (- sz (* 2 dpml)))
(center (+ (* -0.5 sxy) dpml) 0 0)
(eig-match-freq? true)
(eig-parity ODD-Y)
(eig-kpoint (vector3 1.5 0 0))
(eig-resolution 32))))
```

We can exploit the mirror symmetry in the structure and the sources to reduce the computation size by a factor of 2.` ````
(set! symmetries (list (make mirror-sym (direction Y) (phase -1))))
```

We define the near-field surface to span the entire non-PML region above the device, adjacent to the PML in the Z direction.
` ````
(define nearfield
(add-near2far fcen 0 1
(make near2far-region (center 0 0 (- (* 0.5 sz) dpml)) (size (- sxy (* 2 dpml)) (- sxy (* 2 dpml)) 0))))
```

The fields are time stepped until sufficiently decayed.
` ````
(run-sources+ (stop-when-fields-decayed 50 Ey (vector3 0 0 0) 1e-6))
```

Finally, we compute the far fields at 1.55 μm at a series of points on a semicircle within the XZ plane using MEEP's near-to-far-field transformation feature. The radius of the semicircle is 1000 wavelengths which is sufficiently large to obtain the far fields. Each far-field point corresponds to an angle, at equally-spaced intervals, from the X axis. The number of far-field points is a parameter with default of 100. The output consists of the six field components (E` ````
; far-field radius is 1000 wavelengths from the device center
(define-param r (* 1000 (/ fcen)))
; number of far-field points to compute on the semicircle in XZ
(define-param npts 100)
; print the far-field data for each field component at each point on the semicircle
(map (lambda (n)
(let ((ff (get-farfield nearfield (vector3 (* r (cos (* pi (/ n npts)))) 0 (* r (sin (* pi (/ n npts))))))))
(print "farfield:, " (number->string n) ", " (number->string (* pi (/ n npts))))
(map (lambda (m)
(print ", " (number->string (list-ref ff m))))
(arith-sequence 0 1 6))
(print "\n")))
(arith-sequence 0 1 npts))
```

With this simulation script, we can now create an objective function used in the nonlinear optimization for device design. In this example, the objective function takes two parameters as input, the grating periodicity and length (as a fraction of the periodicity), and returns the fraction of the far-field intensity concentrated within the angular cone spanning 70` ````
from numpy import *
from string import *
from subprocess import call
def compute_fieldintensity(p, grad):
a=p[0];
dfrac=p[1];
d=a*dfrac;
np=8;
rot_theta=20;
exec_str = "mpirun -np %d meep-mpi a=%0.2f d=%0.2f rot-theta=%d bragg_outcoupler.ctl > ...
...bragg-optimize-a%0.2f-d%0.2f.out" % (np,a,d,rot_theta,a,d);
call(exec_str, shell="True");
grep_str = "grep farfield: bragg-optimize-a%0.2f-d%0.2f.out |cut -d , -f3- > bragg-optimize-a%0.2f-d%0.2f.dat" % (a,d,a,d);
call(grep_str, shell="True");
mydata = genfromtxt("bragg-optimize-a%0.2f-d%0.2f.dat" % (a,d), delimiter=",", dtype='str');
mydata = char.replace(mydata,'i','j').astype(complex128);
Ex=mydata[:,1]; Ey=mydata[:,2]; Ez=mydata[:,3];
Hx=mydata[:,4]; Hy=mydata[:,5]; Hz=mydata[:,6];
Ex=conj(Ex); Ey=conj(Ey); Ez=conj(Ez);
Px=real(multiply(Ey,Hz)-multiply(Ez,Hy));
Py=real(multiply(Ez,Hx)-multiply(Ex,Hz));
Pz=real(multiply(Ex,Hy)-multiply(Ey,Hx));
Pr=sqrt(square(Px)+square(Py));
Pnorm = Pr/max(Pr);
ang_min = 70;
ang_max = 80;
ang_min = ang_min*pi/180;
ang_max = ang_max*pi/180;
idx = where((mydata[:,0] > ang_min) & (mydata[:,0] < ang_max));
val = sum(Pnorm[idx])/sum(Pnorm);
print "intensity:, %0.2f, %0.2f, %0.6f" % (a,d,val);
return val;
```

Finally, we set up the nonlinear optimization using NLopt by defining the parameter constraints, optimization algorithm, and termination criteria. The initial parameters are chosen randomly. We are using a gradient-free approach since the design involves just two parameters. For designs involving a large number of parameters, a gradient-based approach using adjoint methods would be more efficient. This would enable the use of conjugate gradient methods. The optimization is run multiple times to explore various local optima. A benefit of using NLopt is that we can try several different optimization algorithms by changing just one line in the script. Each solve of the objection function takes approximately 15 minutes on a system with 2.8 GHz AMD Opteron processors.
` ````
import nlopt
import random
from compute_fieldintensity import *
# lattice parameter (um)
a_min = 0.50;
a_max = 3.00;
# waveguide width (fraction of lattice parameter)
dfrac_min = 0.2;
dfrac_max = 0.8;
opt = nlopt.opt(nlopt.LN_BOBYQA, 2)
opt.set_max_objective(compute_fieldintensity)
opt.set_lower_bounds([ a_min, dfrac_min ])
opt.set_upper_bounds([ a_max, dfrac_max ])
opt.set_ftol_abs(0.005)
opt.set_xtol_abs(0.02)
opt.set_initial_step(0.04)
opt.max_eval = 50
# random initial parameters
a_0 = a_min + (a_max-a_min)*random.random();
dfrac_0 = dfrac_min + (dfrac_max-dfrac_min)*random.random();
x = opt.optimize([a_0, dfrac_0])
maxf = opt.last_optimum_value()
print "optimum at a=%0.2f um, d=%0.2f um" % (x[0],x[0]*x[1])
print "maximum value = ", maxf
print "result code = ", opt.last_optimize_result()
```

Here are results from the nonlinear optimization for three different local optima found in the search space. Each run took approximately 10 iterations to converge. The polar plot of the far-field intensity is shown for each design. Also shown is the concentration in the angular cone (objective function). The middle inset has just a single lobe in its radiation pattern and thus the highest concentration.
A schematic of the waveguide unit cell is shown in the figure below. The lattice periodicity (a) is 0.43 μm and the waveguide width (w) and height (h) are 0.50 and 0.22 μm. The hole radius 0.28a which is 0.12 μm. Given the 1d periodicity, we compute the dispersion relation within the irreducible Brillouin zone which spans axial wavevectors along the X direction from 0 to π/a. This is shown in the figure below. There is a bandgap, a region in which there are no guided modes, in the wavelength range of 1.30-1.70 μm. The light line of air is also shown.

` ````
(set-param! resolution 20) ; pixels/a
(define-param a 0.43) ; units of um
(define-param r 0.12) ; units of um
(define-param h 0.22) ; units of um
(define-param w 0.50) ; units of um
(set! r (/ r a)) ; units of "a"
(set! h (/ h a)) ; units of "a"
(set! w (/ w a)) ; units of "a"
(define-param nSi 3.5)
(define Si (make medium (index nSi)))
(set! geometry-lattice (make lattice (size 1 4 4)))
(set! geometry (list (make block (center 0 0 0) (size infinity w h) (material Si))
(make cylinder (center 0 0 0) (radius r) (height infinity) (material air))))
(set! k-points (list (vector3 0 0 0)
(vector3 0.5 0 0)))
(define-param num-kpoints 20)
(set! k-points (interpolate num-kpoints k-points))
(set-param! num-bands 5)
(run-yodd-zeven)
```

We run the simulation script from the shell terminal, pipe the results to a file, and then `grep`

the relevant contents into a separate file for plotting. This takes a few seconds on a machine with a single 2.8 GHz AMD Opteron processor.
` ````
mpb nanobeam-modes.ctl |tee modes.out
grep zevenyoddfreqs: modes.out |cut -d , -f3,7- |sed 1d > modes.dat
```

Finally, we plot the results using Python.
` ````
from numpy import *
import matplotlib.pyplot as plt
f = genfromtxt("modes.dat", delimiter=",")
plt.plot(f[:,0],f[:,1:],'b-',f[:,0],f[:,0],'k-');
plt.xlabel("wavevector k_x (units of 2\pi/a)");
plt.ylabel("frequency (units of 2\pix70^{14} Hz)");
plt.axis([0,0.5,0,0.4]);
plt.show()
```

`a-start`

to ending `a-end`

with radii r=0.28a of several holes on either side of a defect as shown in the schematic below. This cavity design is based on the Applied Physics Letters reference from above. We use MEEP to design the cavity structure to maximize the quality factor (Q) of the fundamental resonant mode. In this example, there is only one design parameter: the cavity length (`s-cav`

).
` ````
(set-param! resolution 40) ; pixels/unit length (1 um)
(define-param a-start 0.43) ; starting periodicity
(define-param a-end 0.33) ; ending periodicity
(define-param s-cav 0.146) ; cavity length
(define-param r 0.28) ; hole radius (units of "a")
(define-param h 0.22) ; waveguide height
(define-param w 0.50) ; waveguide width
(define-param dair 1.00) ; air padding
(define-param dpml 1.00) ; PML thickness
(define-param Ndef 3) ; number of linearly-varying unit cells, excluding starting and ending cells
(define a-taper (interpolate Ndef (list a-start a-end)))
(define a-last (list-ref a-taper (+ Ndef 1)))
(define dgap (- a-last (* 2 r a-last)))
(define-param Nwvg 8) ; number of waveguide unit cells on either side of cavity
(define sx (+ (* 2 (+ (* Nwvg a-start) (fold-left + 0 a-taper))) (- dgap) s-cav))
(define sy (+ dpml dair w dair dpml))
(define sz (+ dpml dair h dair dpml))
(set! geometry-lattice (make lattice (size sx sy sz)))
(set! pml-layers (list (make pml (thickness dpml))))
(define-param nSi 3.5)
(define Si (make medium (index nSi)))
(set! geometry (list (make block (center 0 0 0) (size infinity w h) (material Si))))
(define holes (list '()))
(set! holes (append holes
(map (lambda (mm) (list
(make cylinder (center (+ (* -0.5 sx) (* 0.5 a-start) (* mm a-start)) 0 0) (radius (* r a-start))
(height infinity) (material air))
(make cylinder (center (- (* 0.5 sx) (* 0.5 a-start) (* mm a-start)) 0 0) (radius (* r a-start))
(height infinity) (material air))))
(arith-sequence 0 1 Nwvg))))
(set! holes (apply append holes))
(set! geometry (append geometry holes))
(set! holes (list '()))
(set! holes (append holes
(map (lambda (mm) (list
(make cylinder (center (+ (* -0.5 sx) (* Nwvg a-start)
(if (> mm 0) (fold-left + 0 (map (lambda (nn) (list-ref a-taper nn)) (arith-sequence 0 1 mm))) 0) (* 0.5 (list-ref a-taper mm))) 0 0)
(radius (* r (list-ref a-taper mm))) (height infinity) (material air))
(make cylinder (center (- (* 0.5 sx) (* Nwvg a-start)
(if (> mm 0) (fold-left + 0 (map (lambda (nn) (list-ref a-taper nn)) (arith-sequence 0 1 mm))) 0) (* 0.5 (list-ref a-taper mm))) 0 0)
(radius (* r (list-ref a-taper mm))) (height infinity) (material air))))
(arith-sequence 0 1 (+ Ndef 2)))))
(set! holes (apply append holes))
(set! geometry (append geometry holes))
(define-param lambda-min 1.46) ; minimum source wavelength
(define-param lambda-max 1.66) ; maximum source wavelength
(define fmin (/ lambda-max)) ; minimum source frequency
(define fmax (/ lambda-min)) ; maximum source frequency
(define fcen (* 0.5 (+ fmin fmax))) ; source frequency center
(define df (- fmax fmin)) ; source frequency width
(set! sources (list (make source (src (make gaussian-src (frequency fcen) (fwidth df))) (component Ey) (center 0 0 0))))
(set! symmetries (list (make mirror-sym (direction X) (phase +1))
(make mirror-sym (direction Y) (phase -1))
(make mirror-sym (direction Z) (phase +1))))
(define-param src-time 800)
(run-sources+ src-time (after-sources (harminv Ey (vector3 0 0 0) fcen df))
(in-volume (volume (center 0 0 0) (size sx sy 0)) (at-end output-epsilon output-efield-y)))
```

A bash script is used to sweep the cavity length (`s-cav`

) over the range of 0.110-0.180 μm in increments of 0.005 μm. The output is written to a file with the harminv-related information extracted to a different file for plotting. The figures below show a plot of the resonant wavelength and quality factor as a function of the cavity length. The quality factor is maximum for the design with a cavity length of 0.145 μm. These are consistent with the values reported in the reference.
` ````
#!/bin/bash
for scav in `seq 0.110 0.005 0.180`; do
mpirun -np 2 meep-mpi s-cav=${scav} nanobeam.ctl > nanobeam_cavity_length${scav}.out;
grep harminv0: nanobeam_cavity_length${scav}.out |cut -d , -f2,4 |grep -v frequency >> nanobeam_cavity_varylength.dat;
done;
```

The E` ````
h5topng -o cavity_field_profile.png -C nanobeam-eps-000040000.h5 -vZc bluered nanobeam-ey-000040000.h5
```