This section consists of various projects, with a focus on technology applications, which demonstrate the advanced functionality of MEEP, MPB, and SCUFF/BUFF. For an introduction to each software package including descriptions of core features and application program interface (API) as well as tutorials demonstrating basic functionality, refer to their homepages. Each project consists of a simulation script, a shell script to run the simulation, and iPython and Matlab/Octave post-processing routines for visualizing the output.

MEEP Project #1 — Light-Extraction Efficiency of Organic Light Emitting Diodes (OLEDs)

In this example, we use MEEP to compute the light-extraction efficiency of an organic light-emitting diode (OLED). This is based on results published in Applied Physics Letters, 106, 041111 (2015) (pdf). A typical device structure for a bottom-emitting OLED is shown below. The device consists of a stack of four planar layers. The organic (ORG) layers are deposited on an indium tin oxide (ITO) coated glass substrate with an aluminum (Al) cathode layer on top. Electrons are injected into the organic layer from the Al cathode and holes from the ITO anode. These charge carriers form bound states called excitons which spontaneously recombine to emit photons. Light is extracted from the device through the transparent glass substrate. Some of the light, however, remains trapped within the device as (1) waveguide modes in the high-index ORG/ITO layers and (2) surface-plasmon polaritons (SPP) at the Al/ORG interface. These losses significantly reduce the external quantum efficiency (EQE) of OLEDs. We compute the fraction of the total power in each of these three device components for broadband emission from a white source spanning 400 to 800 nm. The results can be obtained using just one simulation.

There are three key features involved in developing an accurate model. (1) Material Properties. The complex refractive index over the entire broadband spectrum must be imported for each material. This requires fitting the material data to a sum of Drude-Lorentzian susceptibility terms. In this example, we treat the glass, ITO, and organic as lossless since their absorption coefficient is small. The refractive index of Al can be obtained from Applied Optics, 37, pp. 5271-83 (1998). (2) Recombining Excitons as Light Source. An ensemble of spontaneously-recombining excitons produces incoherent emission. This can be modeled using a collection of point-dipole sources with random phase positioned within the organic layer. Given the stochastic nature of the sources, the results must be averaged using Monte-Carlo sampling. The number of samples must be large enough to ensure that the variance of the computed quantities is sufficiently small. (3) Flux Monitors. The total power separated into the three device components is computed using flux monitors. The size and position of these monitors must be chosen correctly to fully capture the relevant fields.

We choose the unit length to be 1 μm. The choice of resolution requires a convergence test to determine a value suitable for a desired level of accuracy. This is particularly relevant in this example given the presence of the lossy-metal aluminum with its nanoscale skin depth as well as the flux monitors. Based on tests, the resolution is set to 100 pixels per μm which is equivalent to a pixel size of 10 nm.
(set-param! resolution 100)          ; pixels/a (a = 1 um)
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. 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)
     (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 in a box configuration to enclose the sources for computing the total power in the device. These monitors must be placed entirely within the organic layer. One flux plane is needed 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)
This simulation script is placed in the file oled-ext-eff.ctl. We use the following Bash script, run_oled.sh, to launch a parallel MEEP simulation using 8 processors.


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 an iPython notebook or Matlab/Octave script. The total power in each device component is the sum of the relevant columns of the output. Here is the Matlab/Octave script with the plots below.
f = load('oled-flux.dat');

% total power emitted from sources
total_power = f(:,2)+f(:,3)+f(:,4)+f(:,5)+f(:,6)+f(:,7);

% total power in glass
glass_power = f(:,8);

% total power in waveguide mode
waveguide_power = f(:,9)+f(:,10)+f(:,11)+f(:,12);

if (length(find(f(:,2:12) < 0)) > 0)
  disp(sprintf("warning: flux is negative"));

% fraction of total power in each of three regions
glass = glass_power./total_power;
waveguide = waveguide_power./total_power;
aluminum = 1-glass-waveguide;

% wavelengths (um)
lambdas = 1./f(:,1);

if (sum(aluminum < 0) > 0)
  disp(sprintf("warning: aluminum absorption is negative"));

if (sum(aluminum > 1) > 0)
  disp(sprintf("warning: aluminum absorption is larger than 1"));

lambdas_linear = linspace(0.4,0.8,100).';
glass_linear = interp1(lambdas,glass,lambdas_linear,"spline","extrap");
waveguide_linear = interp1(lambdas,waveguide,lambdas_linear,"spline","extrap");
aluminum_linear = interp1(lambdas,aluminum,lambdas_linear,"spline","extrap");

xlabel("wavelength (um)"); ylabel("fraction of total power");
legend("glass","aluminum","organic + ITO");
axis([0.4 0.8 0 1]);

disp("power in each region averaged over all wavelengths");
disp(sprintf("glass:     %0.6f, %0.6f",mean(glass_linear),std(glass_linear,1)));
disp(sprintf("aluminum:   %0.6f, %0.6f",mean(aluminum_linear),std(aluminum_linear,1)));
disp(sprintf("organic + ITO: %0.6f, %0.6f",mean(waveguide_linear),std(waveguide_linear,1)));			    
There are significant losses from the SPPs for the perpendicular dipoles as shown in the figure on the right. This is because this dipole orientation couples readily into these modes due to matching boundary conditions of the electric fields. Such losses are smaller for the parallel dipoles. The total power extracted from the device is computing using ⅔P+⅓P≈ 0.20. This value is consistent with experimentally-measured results. As described in the reference paper, we can apply a nanoscale surface texture to the cathode in order to recover the light lost to the SPPs, which dominate the total losses, and also to enhance the rate of spontaneous emission via the Purcell effect.