Meep Tutorial/Optical forces

From AbInitio

(Difference between revisions)
Jump to: navigation, search
Revision as of 08:44, 1 October 2012 (edit)
Ardavan (Talk | contribs)

← Previous diff
Revision as of 08:57, 1 October 2012 (edit)
Ardavan (Talk | contribs)

Next diff →
Line 5: Line 5:
:<math>F=-\frac{1}{\omega}\frac{d\omega}{ds}\Bigg\vert_\vec{k}U,</math> :<math>F=-\frac{1}{\omega}\frac{d\omega}{ds}\Bigg\vert_\vec{k}U,</math>
-where <math>\omega</math> is the eigenmode frequency of the coupled-waveguide system, <math>s</math> is the separation distance between the parallel waveguides, <math>k</math> is the conserved wave vector and <math>U</math> is the total energy of the electromagnetic fields. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. This expression has been shown to be mathematically equivalent to the Maxwell stress tensor method and in this tutorial we will verify this result for the parallel-waveguide example. It is convenient to normalize the force so as to eliminate units altogether and since we know that the total power transmitted through the waveguide is <math>P=v_gU/L</math> (<math>v_g</math> is the group velocity, <math>L</math> is the waveguide length and <math>U</math> is defined as before), we focus instead on the dimensionless quantity force per unit length and power <math>(F/L)(ac/P)</math> (<math>a</math> is an arbitrary unit length and <math>c</math> is the speed of light).+where <math>\omega</math> is the eigenmode frequency of the coupled-waveguide system, <math>s</math> is the separation distance between the parallel waveguides, <math>k</math> is the conserved wave vector and <math>U</math> is the total energy of the electromagnetic fields. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. This expression has been shown to be mathematically equivalent to the Maxwell stress tensor method and in this tutorial we will verify this result for the parallel-waveguide example. It is convenient to normalize the force so as to eliminate units altogether and since we know that the total power transmitted through the waveguide is <math>P=v_gU/L</math> (<math>v_g</math> is the group velocity, <math>L</math> is the waveguide length and <math>U</math> is defined as before), we focus instead on the dimensionless quantity: force per unit length and power <math>(F/L)(ac/P)</math> (<math>a</math> is an arbitrary unit length and <math>c</math> is the speed of light).
-We can therefore compute the optical gradient force in two ways, one indirect and the other direct: 1) in the frequency-domain planewave-solver MPB, we compute the eigenmode frequency and group velocity for a mode with a fixed propagating wavevector for a number of different separation distances and then use a finite-difference scheme to evaluate the above expression and 2) in Meep, we compute the force directly using the MST as well as the power transmitted through the waveguide only at the frequency corresponding to the guided mode. In this particular example, we consider only the fundamental <code>y-odd</code> mode, which happens to clearly show the bi-directional feature of the force.+We can therefore compute the optical gradient force in two ways, one indirect and the other direct: 1) using the frequency-domain planewave-solver MPB, we compute the eigenmode frequency (and group velocity) for a mode with a fixed propagating wavevector at a number of different separation distances and then use a finite-difference scheme to evaluate the above expression and 2) in Meep, we compute the force directly using the MST as well as the power transmitted through the waveguide only at the frequency corresponding to the guided mode. In this particular example, we consider just the fundamental <code>y-odd</code> mode, which happens to clearly show the bi-directional features of the force.
First, let us set up the two-dimensional computational cell (although since we are interested in a propagating mode of the structure which requires an axial wavevector, technically we will have a three-dimensional simulation) First, let us set up the two-dimensional computational cell (although since we are interested in a propagating mode of the structure which requires an axial wavevector, technically we will have a three-dimensional simulation)
Line 14: Line 14:
(define-param nSi 3.45) (define-param nSi 3.45)
(define Si (make medium (index nSi))) (define Si (make medium (index nSi)))
- (define-param dpml 5.0)+ (define-param dpml 1.0)
- (define-param sx 20)+ (define-param sx 5)
- (define-param sy 5)+ (define-param sy 3)
(set! geometry-lattice (set! geometry-lattice
(make lattice (size (+ sx (* 2 dpml)) (+ sy (* 2 dpml)) no-size))) (make lattice (size (+ sx (* 2 dpml)) (+ sy (* 2 dpml)) no-size)))
(set! pml-layers (list (make pml (thickness dpml)))) (set! pml-layers (list (make pml (thickness dpml))))
(define-param sw 1.0) ; waveguide width (define-param sw 1.0) ; waveguide width
- (define-param d 1.0) ; waveguide separation distance+ (define-param s 1.0) ; waveguide separation distance
(set! geometry (list (set! geometry (list
- (make block (center (* -0.5 (+ d sw)) 0)+ (make block (center (* -0.5 (+ s sw)) 0)
(size sw sw infinity) (material Si)) (size sw sw infinity) (material Si))
- (make block (center (* 0.5 (+ d sw)) 0)+ (make block (center (* 0.5 (+ s sw)) 0)
(size sw sw infinity) (material Si)))) (size sw sw infinity) (material Si))))
Line 36: Line 36:
(make mirror-sym (direction Y) (phase (if yodd? -1 +1))))) (make mirror-sym (direction Y) (phase (if yodd? -1 +1)))))
-Next, we set the Bloch-periodic boundary condition in order to excite the specific guided mode of the waveguide system+Next, we set the Bloch-periodic boundary condition in order to excite the specific guided mode of the waveguide system (corresponding to a wavevector of <math>\pi/a</math>)
(define-param beta 0.5) (define-param beta 0.5)
(set! k-point (vector3 0 0 beta)) (set! k-point (vector3 0 0 beta))
-Since we do not know a-priori what the eigenmode frequency will be for the given separation distance, we first use a point-dipole source positioned in the middle of the waveguides with a Gaussian profile and then use <code>harminv</code> to determine the resonant mode frequency.+Since we do not know a-priori what the eigenmode frequency will be at any given separation distance, we first use a point-dipole source positioned in the middle of the waveguides with a Gaussian profile and then use <code>harminv</code> to determine the resonant mode frequency.
(define-param fcen 0.2) (define-param fcen 0.2)
Line 47: Line 47:
(set! sources (list (set! sources (list
(make source (src (make gaussian-src (frequency fcen) (fwidth df))) (make source (src (make gaussian-src (frequency fcen) (fwidth df)))
- (component Ey) (center (* -0.5 (+ d sw)) 0))+ (component Ey) (center (* -0.5 (+ s sw)) 0))
(make source (src (make gaussian-src (frequency fcen) (fwidth df))) (make source (src (make gaussian-src (frequency fcen) (fwidth df)))
- (component Ey) (center (* 0.5 (+ d sw)) 0) (amplitude (if xodd? -1.0 1.0)))))+ (component Ey) (center (* 0.5 (+ s sw)) 0) (amplitude (if xodd? -1.0 1.0)))))
(define-param runtime 300) (define-param runtime 300)
(run-sources+ runtime (run-sources+ runtime
- (after-sources (harminv Ey (vector3 (* 0.5 (+ d sw)) 0) fcen df)))+ (after-sources (harminv Ey (vector3 (* 0.5 (+ s sw)) 0) fcen df)))
(define f (harminv-freq-re (car harminv-results))) (define f (harminv-freq-re (car harminv-results)))
- (print "freq:, " d ", " f "\n")+ (print "freq:, " s ", " f "\n")
-Once we have determined what the eigenmode frequency is, we can then use this value to precisely excite the mode of interest using the the <code>eigenmode-source</code> feature and also compute the force and the transmitted power only at this value.+Once we have determined what the eigenmode frequency is, we can then use this value to accurately excite the mode of interest using the <code>eigenmode-source</code> feature and also compute the force and the transmitted power only at this value.
(reset-meep) (reset-meep)
Line 64: Line 64:
(component Ey) (component Ey)
(size sw sw 0) (size sw sw 0)
- (center (* -0.5 (+ d sw)) 0)+ (center (* -0.5 (+ s sw)) 0)
(eig-kpoint k-point) (eig-kpoint k-point)
(eig-match-freq? true) (eig-match-freq? true)
Line 72: Line 72:
(component Ey) (component Ey)
(size sw sw 0) (size sw sw 0)
- (center (* 0.5 (+ d sw)) 0)+ (center (* 0.5 (+ s sw)) 0)
(eig-kpoint k-point) (eig-kpoint k-point)
(eig-match-freq? true) (eig-match-freq? true)
Line 79: Line 79:
(define wvg-pwr (add-flux f 0 1 (define wvg-pwr (add-flux f 0 1
(make flux-region (direction Z) (center 0 0) (make flux-region (direction Z) (center 0 0)
- (size (* 1.2 (+ (* 2 sw) d)) (* 1.2 sw) 0))))+ (size (* 1.2 (+ (* 2 sw) s)) (* 1.2 sw) 0))))
(define-param dpad 0.1) (define-param dpad 0.1)
(define wvg-force (add-force f 0 1 (define wvg-force (add-force f 0 1
(make force-region (direction X) (weight +1.0) (make force-region (direction X) (weight +1.0)
- (center (- (* 0.5 d) dpad) 0) (size 0 sy))+ (center (- (* 0.5 s) dpad) 0) (size 0 sy))
(make force-region (direction X) (weight -1.0) (make force-region (direction X) (weight -1.0)
- (center (+ (* 0.5 d) sw dpad) 0) (size 0 sy))))+ (center (+ (* 0.5 s) sw dpad) 0) (size 0 sy))))
(run-sources+ runtime) (run-sources+ runtime)
(display-fluxes wvg-pwr) (display-fluxes wvg-pwr)
(display-forces wvg-force) (display-forces wvg-force)
-Note that we have defined the flux plane to contain an area spanning both waveguides rather than two separate flux planes, one for each waveguide. This is because in the limit of small separation, two separate flux planes will overlap thus complicating the analysis whereas the total power through a single, all-encompassing flux plane need only be halved in order to determine the value for just one of the two waveguides. Also note that instead of defining a closed, four-sided box surrounding the waveguide when computing the MST, we chose instead to compute the MST along two vertically-oriented lines, one on each side of the waveguide since by symmetry we need not consider the force in the <code>y</code> direction. +Note that we have defined the flux plane to span an area slightly larger than both waveguides rather than two separate flux planes, one for each waveguide. This is because in the limit of small separation, two separate flux planes will overlap thus complicating the analysis whereas the total power through a single, all-encompassing flux plane need only be halved in order to determine the value for just one of the two waveguides. Also note that instead of defining a closed, four-sided box surrounding the waveguide when computing the MST, we chose instead to compute the MST along two vertically-oriented lines, one on each side of the waveguide since by symmetry we need not consider the force in the <code>y</code> direction.
[[Image:Waveguide_forces.png|center|Normalized force per unit length and input power versus waveguide separation of two parallel waveguides.]] [[Image:Waveguide_forces.png|center|Normalized force per unit length and input power versus waveguide separation of two parallel waveguides.]]

Revision as of 08:57, 1 October 2012

Here we will demonstrate Meep's ability to compute classical forces using the Maxwell stress tensor (MST) as well as the eigenmode source feature which integrates the frequency-domain functionality of our other open-source tool MPB. Our example consists of two identical dielectric waveguides (made of non-lossy silicon with a square cross section) supporting the same mode while separated in space. Due to the symmetry and orientation of the waveguide profile, the two modes can be chosen to be either symmetric or anti-symmetric with respect to a mirror plane between them. As the two waveguides are brought closer and closer together, the individual guided modes interact and give rise to an optical gradient force that is transverse to the waveguide axis (this is to be contrasted with radiation pressure that involves momentum exchange between photons and is longitudinal in nature). An interesting phenomena that occurs is that the sign of the force can be tuned to be either attractive or repulsive depending on just the relative phase of the two modes which we will demonstrate in this tutorial.

The optical gradient force arising from the evanescent coupling of the modes of two adjacent structures can, in addition to the MST, be computed using the following analytic expression:

F=-\frac{1}{\omega}\frac{d\omega}{ds}\Bigg\vert_\vec{k}U,

where ω is the eigenmode frequency of the coupled-waveguide system, s is the separation distance between the parallel waveguides, k is the conserved wave vector and U is the total energy of the electromagnetic fields. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. This expression has been shown to be mathematically equivalent to the Maxwell stress tensor method and in this tutorial we will verify this result for the parallel-waveguide example. It is convenient to normalize the force so as to eliminate units altogether and since we know that the total power transmitted through the waveguide is P = vgU / L (vg is the group velocity, L is the waveguide length and U is defined as before), we focus instead on the dimensionless quantity: force per unit length and power (F / L)(ac / P) (a is an arbitrary unit length and c is the speed of light).

We can therefore compute the optical gradient force in two ways, one indirect and the other direct: 1) using the frequency-domain planewave-solver MPB, we compute the eigenmode frequency (and group velocity) for a mode with a fixed propagating wavevector at a number of different separation distances and then use a finite-difference scheme to evaluate the above expression and 2) in Meep, we compute the force directly using the MST as well as the power transmitted through the waveguide only at the frequency corresponding to the guided mode. In this particular example, we consider just the fundamental y-odd mode, which happens to clearly show the bi-directional features of the force.

First, let us set up the two-dimensional computational cell (although since we are interested in a propagating mode of the structure which requires an axial wavevector, technically we will have a three-dimensional simulation)

 (set-param! resolution 20)
 (define-param nSi 3.45)
 (define Si (make medium (index nSi)))
 (define-param dpml 1.0)
 (define-param sx 5)
 (define-param sy 3)
 (set! geometry-lattice 
     (make lattice (size (+ sx (* 2 dpml)) (+ sy (* 2 dpml)) no-size)))
 (set! pml-layers (list (make pml (thickness dpml))))
 (define-param sw 1.0) ; waveguide width
 (define-param s 1.0)  ; waveguide separation distance
 (set! geometry (list 
       (make block (center (* -0.5 (+ s sw)) 0)
             (size sw sw infinity) (material Si))
       (make block (center (*  0.5 (+ s sw)) 0)
             (size sw sw infinity) (material Si))))

There are two mirror symmetries that we can exploit to reduce the simulation size by a factor of four.

 (define-param xodd? true)
 (define-param yodd? true)
 (set! symmetries (list 
        (make mirror-sym (direction X) (phase (if xodd? -1 +1)))
        (make mirror-sym (direction Y) (phase (if yodd? -1 +1)))))

Next, we set the Bloch-periodic boundary condition in order to excite the specific guided mode of the waveguide system (corresponding to a wavevector of π / a)

 (define-param beta 0.5)
 (set! k-point (vector3 0 0 beta))

Since we do not know a-priori what the eigenmode frequency will be at any given separation distance, we first use a point-dipole source positioned in the middle of the waveguides with a Gaussian profile and then use harminv to determine the resonant mode frequency.

 (define-param fcen 0.2)
 (define-param df 0.1)
 (set! sources (list 
         (make source (src (make gaussian-src (frequency fcen) (fwidth df))) 
                (component Ey) (center (* -0.5 (+ s sw)) 0))
         (make source (src (make gaussian-src (frequency fcen) (fwidth df))) 
                (component Ey) (center (* 0.5 (+ s sw)) 0) (amplitude (if xodd? -1.0 1.0)))))
 (define-param runtime 300)
 (run-sources+ runtime 
     (after-sources (harminv Ey (vector3 (* 0.5 (+ s sw)) 0) fcen df)))
 (define f (harminv-freq-re (car harminv-results)))
 (print "freq:, " s ", " f "\n")

Once we have determined what the eigenmode frequency is, we can then use this value to accurately excite the mode of interest using the eigenmode-source feature and also compute the force and the transmitted power only at this value.

 (reset-meep)
 (change-sources! (list
    (make eigenmode-source
      (src (make gaussian-src (frequency f) (fwidth df)))
      (component Ey)
      (size sw sw 0)
      (center (* -0.5 (+ s sw)) 0)
      (eig-kpoint k-point)
      (eig-match-freq? true)
      (eig-parity (if yodd? ODD-Y EVEN-Y)))
    (make eigenmode-source
      (src (make gaussian-src (frequency f) (fwidth df)))
      (component Ey)
      (size sw sw 0)
      (center (* 0.5 (+ s sw)) 0)
      (eig-kpoint k-point)
      (eig-match-freq? true)
      (eig-parity (if yodd? ODD-Y EVEN-Y))
      (amplitude (if xodd? -1.0 1.0)))))
 (define wvg-pwr (add-flux f 0 1
    (make flux-region (direction Z) (center 0 0) 
 	   (size (* 1.2 (+ (* 2 sw) s)) (* 1.2 sw) 0))))
 (define-param dpad 0.1)
 (define wvg-force (add-force f 0 1
   (make force-region (direction X) (weight +1.0) 
         (center (- (* 0.5 s) dpad) 0) (size 0 sy))
   (make force-region (direction X) (weight -1.0) 
         (center (+ (* 0.5 s) sw dpad) 0) (size 0 sy))))
 (run-sources+ runtime)
 (display-fluxes wvg-pwr)
 (display-forces wvg-force)

Note that we have defined the flux plane to span an area slightly larger than both waveguides rather than two separate flux planes, one for each waveguide. This is because in the limit of small separation, two separate flux planes will overlap thus complicating the analysis whereas the total power through a single, all-encompassing flux plane need only be halved in order to determine the value for just one of the two waveguides. Also note that instead of defining a closed, four-sided box surrounding the waveguide when computing the MST, we chose instead to compute the MST along two vertically-oriented lines, one on each side of the waveguide since by symmetry we need not consider the force in the y direction.

Normalized force per unit length and input power versus waveguide separation of two parallel waveguides.
Personal tools