A 2d MEEP simulation of the transmission power around a 90-degree Si waveguide bend. The excitation source is TM polarized (i.e. E-field is out of plane).
Here we would like to know exactly how much power makes it around a 90 degrees Si waveguide bend, how much is reflected, and how much is radiated away. How can we do this?
Software Used
MEEP
Meep (or MEEP) is a finite-difference time-domain (FDTD) simulation software package developed at MIT to model electromagnetic systems, along with the MPB eigenmode package which is also available on Kogence. MEEP stands for MIT Electromagnetic Equation Propagation.
For more details see MEEP .
Disclaimer
The material in the discussion below was taken from the MEEP tutorial in [1].
Input File
We'll tell MEEP to keep track of the fields and their Fourier transforms in a certain region, and from this compute the flux of electromagnetic energy as a function of ω. Moreover, we'll get an entire spectrum of the transmission in a single run, by Fourier-transforming the response to a short pulse. However, in order to normalize the transmission (to get transmission as a fraction of incident power), we'll have to do two runs, one with and one without a bend.
We first define the cell size as parameters:
(define-param sx 16) ; size of cell in X direction
(define-param sy 32) ; size of cell in Y direction
(set! geometry-lattice (make lattice (size sx sy no-size)))
wheredefine-param
is a libctl feature to define variables that can be overridden from the command line. We could now do meep sx=17 tut-wvg-bend-trans.ctl
to change the X size to 17, without editing the ctl file, for example. We'll also define a couple of parameters to set the width of the waveguide and the "padding" between it and the edge of the computational cell:
(define-param pad 4) ; padding distance between waveguide and cell edge
(define-param w 1) ; width of waveguide
In order to define the waveguide positions, etcetera, we will now have to use arithmetic. For example, the y center of the horizontal waveguide will be given by -0.5 * (sy - w - 2*pad)
. At least, that is what the expression would look like in C; in Scheme, the syntax for 1 + 2 is (+ 1 2)
, and so on, so we will define the vertical and horizontal waveguide centers as:
(define wvg-ycen (* -0.5 (- sy w (* 2 pad)))) ; y center of horiz. wvg
(define wvg-xcen (* 0.5 (- sx w (* 2 pad)))) ; x center of vert. wvg
Now, we have to make the geometry, as before. This time, however, we really want two geometries: the bend, and also a straight waveguide for normalization. We could do this with two separate ctl files, but that is annoying. Instead, we'll define a parameter no-bend?
which is true
for the straight-waveguide case and false
for the bend.
(define-param no-bend? false) ; if true, have straight waveguide, not bend
Now, we define the geometry via two cases, with an if statement—the Scheme syntax is (if predicate? if-true if-false)
.
1 (set! geometry
2 (if no-bend?
3 (list
4 (make block
5 (center 0 wvg-ycen)
6 (size infinity w infinity)
7 (material (make dielectric (epsilon 12)))))
8 (list
9 (make block
10 (center (* -0.5 pad) wvg-ycen)
11 (size (- sx pad) w infinity)
12 (material (make dielectric (epsilon 12))))
13 (make block
14 (center wvg-xcen (* 0.5 pad))
15 (size w (- sy pad) infinity)
16 (material (make dielectric (epsilon 12)))))))
Thus, if no-bend?
is true
we make a single block for a straight waveguide, and otherwise we make two blocks for a bent waveguide.
The source is now a gaussian-src
instead of a continuous-src
, parameterized by a center frequency and a frequency width (the width of the Gaussian spectrum), which we'll define via define-param
as usual.
1 (define-param fcen 0.15) ; pulse center frequency
2 (define-param df 0.1) ; pulse width (in frequency)
3 (set! sources (list
4 (make source
5 (src (make gaussian-src (frequency fcen) (fwidth df)))
6 (component Ez)
7 (center (+ 1 (* -0.5 sx)) wvg-ycen)
8 (size 0 w))))
Notice how we're using our parameters like wvg-ycen
and w
: if we change the dimensions, everything will now shift automatically. The boundary conditions and resolution are set as before, except that now we'll use set-param!
so that we can override the resolution from the command line.:
(set! pml-layers (list (make pml (thickness 1.0))))
(set-param! resolution 10)
Finally, we have to specify where we want MEEP to compute the flux spectra, and at what frequencies. (This must be done after specifying the geometry, sources, resolution, etcetera, because all of the field parameters are initialized when flux planes are created.)
1 (define-param nfreq 100) ; number of frequencies at which to compute flux
2 (define trans ; transmitted flux
3 (add-flux fcen df nfreq
4 (if no-bend?
5 (make flux-region
6 (center (- (/ sx 2) 1.5) wvg-ycen) (size 0 (* w 2)))
7 (make flux-region
8 (center wvg-xcen (- (/ sy 2) 1.5)) (size (* w 2) 0)))))
9 (define refl ; reflected flux
10 (add-flux fcen df nfreq
11 (make flux-region
12 (center (+ (* -0.5 sx) 1.5) wvg-ycen) (size 0 (* w 2)))))
We compute the fluxes through a line segment twice the width of the waveguide, located at the beginning or end of the waveguide. (Note that the flux lines are separated by 1 from the boundary of the cell, so that they do not lie within the absorbing PML regions.) Again, there are two cases: the transmitted flux is either computed at the right or the bottom of the computational cell, depending on whether the waveguide is straight or bent.
Here, the fluxes will be computed for 100 (nfreq
) frequencies centered on fcen
, from fcen-df/2
to fcen+df/2
. That is, we only compute fluxes for frequencies within our pulse bandwidth. This is important because, to far outside the pulse bandwidth, the spectral power is so low that numerical errors make the computed fluxes useless.
Result
Now computing reflection spectra is a bit tricky because we need to separate the incident and reflected fields. We do this in MEEP by saving the Fourier-transformed fields from the normalization run (no-bend?=true
), and loading them, negated, before the other runs. The latter subtracts the Fourier-transformed incident fields from the Fourier transforms of the scattered fields; logically, we might subtract these after the run, but it turns out to be more convenient to subtract the incident fields first and then accumulate the Fourier transform. All of this is accomplished with two commands, save-flux
(after the normalization run) and load-minus-flux
(before the other runs). We can call them as follows:
(if (not no-bend?) (load-minus-flux "refl-flux" refl))
(run-sources+ 500 (at-beginning output-epsilon))
(if no-bend? (save-flux "refl-flux" refl))
This uses a file called bend-flux-refl-flux.h5
(the ctl file name is used as a prefix) to store/load the Fourier transformed fields in the flux planes. The (run-sources+ 500)
runs the simulation until the Gaussian source has turned off (which is done automatically once it has decayed for a few standard deviations), plus an additional 500 time units.
Why do we keep running after the source has turned off? Because we must give the pulse time to propagate completely across the cell. Moreover, the time required is a bit tricky to predict when you have complex structures, because there might be resonant phenomena that allow the source to bounce around for a long time. Therefore, it is convenient to specify the run time in a different way: instead of using a fixed time, we require that the | E_{z} | ^{2} at the end of the waveguide must have decayed by a given amount (e.g. 1/1000) from its peak value. We can do this via:
1 (run-sources+
2 (stop-when-fields-decayed 50 Ez
3 (if no-bend?
4 (vector3 (- (/ sx 2) 1.5) wvg-ycen)
5 (vector3 wvg-xcen (- (/ sy 2) 1.5)))
6 0.001))
stop-when-fields-decayed
takes four arguments: (stop-when-fields-decayed dT component pt decay-by)
. What it does is, after the sources have turned off, it keeps running for an additional dT
time units every time the given |component|^{2} at the given point has not decayed by at least decay-by
from its peak value for all times within the previous dT
. In this case, dT=50
, the component is E_{z}, the point is at the center of the flux plane at the end of the waveguide, and decay-by=0.001
. So, it keeps running for an additional 50 time units until the square amplitude has decayed by 1/1000 from its peak: this should be sufficient to ensure that the Fourier transforms have converged.
Finally, we have to output the flux values:
(display-fluxes trans refl)
This prints a series of outputs like:
flux1:, 0.1, 7.91772317108475e-7, -3.16449591437196e-7
flux1:, 0.101010101010101, 1.18410865137737e-6, -4.85527604203706e-7
flux1:, 0.102020202020202, 1.77218779386503e-6, -7.37944901819701e-7
flux1:, 0.103030303030303, 2.63090852112034e-6, -1.11118350510327e-6
flux1:, ...
This is comma-delimited data, which can easily be imported into any spreadsheet or plotting program (e.g. GNUPlot): the first column is the frequency, the second is the transmitted power, and the third is the reflected power.
Now, we need to run the simulation twice, once with no-bend?=true
and once with no-bend?=false
(the default):
meep no-bend?=true bend-flux.ctl | tee bend0.out
meep bend-flux.ctl | tee bend.out
(The tee
command is a useful Unix command that saves the output to a file and displays it on the screen, so that we can see what is going on as it runs.) Then, we should pull out the flux1
lines into a separate file to import them into our plotting program:
grep flux1: bend0.out > bend0.dat unix% grep flux1: bend.out > bend.dat
We next plot the results using GNUPlot:
What are we plotting here? The transmission is the transmitted flux (second column of bend.dat
) divided by the incident flux (second column of bend0.dat
), to give us the fraction of power transmitted. The reflection is the reflected flux (third column of bend.dat
) divided by the incident flux (second column of bend0.dat
); we also have to multiply by −1 because all fluxes in Meep are computed in the positive-coordinate direction by default, and we want the flux in the −x direction. Finally, the loss is simply:
1 - transmission - reflection.
We should also check whether our data is converged, by increasing the resolution and cell size and seeing by how much the numbers change. In this case, we'll just try doubling the cell size:
meep sx=32 sy=64 no-bend?=true bend-flux.ctl |tee bend0-big.out
meep sx=32 sy=64 bend-flux.ctl |tee bend-big.out
Again, we must run both simulations in order to get the normalization right. The results are included in the plot above as grey dashed lines—you can see that the numbers have changed slightly for transmission and loss, probably stemming from interference between light radiated directly from the source and light propagating around the waveguide. To be really confident, we should probably run the simulation again with an even bigger cell, but we'll call it enough for this tutorial.
Related Models
Following are some related models available for cloning/copying by anyone:
- Calculate transmission flux around a 90-degree Si waveguide bend in 2d. mohan12345
- Calculate transmission flux around a 90-degree Si waveguide bend in 2d. testing
- Calculate transmission flux around a 90-degree Si waveguide bend in 2d. mohan
- Calculate transmission flux around a 90-degree Si waveguide bend in 2d.
- Calculate transmission flux around a 90-degree Si waveguide bend in 2d. copying windows
- Calculate 2d ring-resonator TM modes using cartesian coordinates
- TE Band structure of a Si PC waveguide with holes
- Field profile in a 90-degree Si waveguide bend using 2d FDTD - MEEP
- Calculate 2d ring-resonator TM modes using cylindrical coordinates
- Calculate 2d ring-resonator TM modes using cylindrical coordinates - my copy1
- Calculate 2d ring-resonator TM modes using cylindrical coordinates - joaquin
Click on the category links at the bottom of this page to navigate to a full list of simulation models in similar subject area or similar computational methodology.
References
[1] http://ab-initio.mit.edu/wiki/index.php/Meep_tutorial#Transmission_spectrum_around_a_waveguide_bend