Meep Tutorial/Optical forces

From AbInitio

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

← Previous diff
Current revision (01:50, 3 May 2017) (edit)
Ardavan (Talk | contribs)

 
Line 1: Line 1:
-Here we will demonstrate Meep's ability to compute classical forces using the [[w: Maxwell stress tensor|Maxwell stress tensor]] (MST) as well as the eigenmode source feature which integrates the frequency-domain functionality of our other open-source tool [[MPB]] (note that this will require Meep 1.2 or greater). Our example consists of two identical dielectric waveguides (made of non-lossy silicon with a square cross section) supporting an identical 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 more and more 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 for this particular geometry 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.+Here we will demonstrate Meep's ability to compute classical forces using the [[w: Maxwell stress tensor|Maxwell stress tensor]] (MST) as well as the eigenmode source feature which integrates our mode-solver package [[MPB]]. Note that this will require Meep 1.2+ and MPB 1.5+. Our example consists of two identical dielectric waveguides made of lossless Silicon with a square cross section supporting identical modes and separated in air. Due to the symmetry and orientation of the waveguides, 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 couple more and more and give rise to an optical gradient force that is <i>transverse</i> to the waveguide axis. This is to be contrasted with [[w: Radiation pressure|radiation pressure]] that involves momentum exchange between photons and is <i>longitudinal</i> in nature. An interesting phenomena that occurs for this system is that the <i>sign</i> of the force can be tuned to be either attractive or repulsive depending on the relative phase of the two modes. We will demonstrate this effect 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: 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:
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 (for more details, see the original [http://math.mit.edu/~stevenj/papers/PovinelliLo05.pdf paper] by Povinelli et al.). 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 [http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-20-18116 result] by Rakich et al. for the parallel-waveguide example. It is convenient to normalize the force so as to eliminate the tricky 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 following 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. For more details, refer to [http://math.mit.edu/~stevenj/papers/PovinelliLo05.pdf Povinelli et al.]. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. This expression has been shown to be mathematically equivalent to the MST (see [http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-20-18116 Rakich et al.]) 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 the tricky units altogether. Since the total power transmitted through the waveguide is <math>P=v_gU/L</math> where <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 force per unit length and power <math>(F/L)(ac/P)</math> where <math>a</math> is an arbitrary unit length and <math>c</math> is the speed of light. This dimensionless quantity enables us to compute both the flux and the force in a single simulation.
-We can therefore compute the optical gradient force in two ways, one indirect and the other direct: 1) using any frequency-domain solver (in our case MPB), we compute the eigenmode frequency (and corresponding 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 both the force using the MST and 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.+We can therefore compute the optical gradient force in two ways, one indirect and the other direct: 1) using any mode solver (in our case MPB), we compute the eigenmode frequency (and corresponding 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 both the force using the MST and 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 defined later on, technically we will have a three-dimensional simulation)+First, let's set up the two-dimensional computational cell (although since we are interested in a propagating mode of the structure which requires an axial wavevector, this will technically be a three-dimensional simulation)
(set-param! resolution 30) (set-param! resolution 30)
Line 20: Line 20:
(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 a 1.0) ; waveguide width
(define-param s 1.0) ; waveguide separation distance (define-param s 1.0) ; waveguide separation distance
(set! geometry (list (set! geometry (list
- (make block (center (* -0.5 (+ s sw)) 0)+ (make block (center (* -0.5 (+ s a)) 0)
- (size sw sw infinity) (material Si))+ (size a a infinity) (material Si))
- (make block (center (* 0.5 (+ s sw)) 0)+ (make block (center (* 0.5 (+ s a)) 0)
- (size sw sw infinity) (material Si))))+ (size a a infinity) (material Si))))
There are two mirror symmetries that we can exploit to reduce the simulation size by a factor of four. There are two mirror symmetries that we can exploit to reduce the simulation size by a factor of four.
Line 40: Line 40:
(set! k-point (vector3 0 0 beta)) (set! k-point (vector3 0 0 beta))
-Since we do not know apriori what the eigenmode frequency will be at any given separation distance, we first excite a spectrum of frequencies using a broadband point-dipole source positioned in the middle of each waveguides and then determine the resonant mode frequency using <code>harminv</code>. Since we are using a Bloch periodic boundary condition in the <code>z</code> direction, the excited mode (if any) will propagate indefinitely in time which is why we stop the simulation at a fixed <code>runtime</code> after the Gaussian sources have turned off (choosing a suitable <code>runtime</code> requires some care since larger values will in turn lead to larger fields which might cause instabilities in the floating-point arithmetic).+Since we do not know apriori what the eigenmode frequency will be at any given separation distance, we first excite a spectrum of frequencies using a broadband volume source in each waveguide and then determine the resonant frequency using <code>harminv</code>. The use of a Bloch periodic boundary condition in the <code>z</code> direction will mean that the excited mode (if any) will propagate indefinitely in time which is why we stop the simulation at a fixed 200 time units after the Gaussian-pulsed sources have turned off.
(define-param fcen 0.22) (define-param fcen 0.22)
Line 46: Line 46:
(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 (+ s sw)) 0))+ (component Ey) (center (* -0.5 (+ s a)) 0) (size a a 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 (+ s sw)) 0) (amplitude (if xodd? -1.0 1.0)))))+ (component Ey) (center (* 0.5 (+ s a)) 0) (size a a 0)
- (define-param runtime 200)+ (amplitude (if xodd? -1.0 1.0)))))
- (run-sources+ runtime + (run-sources+ 200
- (after-sources (harminv Ey (vector3 (* 0.5 (+ s sw)) 0) fcen df)))+ (after-sources (harminv Ey (vector3 (* 0.5 (+ s a)) 0) fcen df)))
(define f (harminv-freq-re (car harminv-results))) (define f (harminv-freq-re (car harminv-results)))
(print "freq:, " s ", " f "\n") (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 <code>eigenmode-source</code> feature and also compute the force and the transmitted power only at this value. The <code>eigenmode-mode</code> feature invokes [[MPB]] via a library routine in order to compute the relevant mode of interest and subsequently imports the steady-state field profile into the Meep simulation for use as the initial amplitude of the source. This then enables an efficient excitation of the relevant mode of interest in the time domain to a much higher degree of accuracy than would otherwise be possible had we simply used a point-dipole source.+Once we have determined what the eigenmode frequency is, we then use this value to accurately excite the mode of interest using the <code>eigenmode-source</code> feature and also compute the force on the waveguide and the power in the guided mode only at this value. The <code>eigenmode-mode</code> feature invokes [[MPB]] via a library routine in order to compute the relevant mode of interest and subsequently imports the (steady-state) field profile into the Meep simulation for use as the initial amplitude of the source. This enables an efficient excitation of the relevant eigenmode of interest in the time domain to a much higher degree of accuracy than would otherwise be possible had we simply used a point-dipole source. For more details, refer to [http://arxiv.org/abs/1301.5366 chapter 4 of our book].
(reset-meep) (reset-meep)
Line 62: Line 62:
(src (make gaussian-src (frequency f) (fwidth df))) (src (make gaussian-src (frequency f) (fwidth df)))
(component Ey) (component Ey)
- (size sw sw 0)+ (size a a 0)
- (center (* -0.5 (+ s sw)) 0)+ (center (* -0.5 (+ s a)) 0)
(eig-kpoint k-point) (eig-kpoint k-point)
(eig-match-freq? true) (eig-match-freq? true)
Line 70: Line 70:
(src (make gaussian-src (frequency f) (fwidth df))) (src (make gaussian-src (frequency f) (fwidth df)))
(component Ey) (component Ey)
- (size sw sw 0)+ (size a a 0)
- (center (* 0.5 (+ s sw)) 0)+ (center (* 0.5 (+ s a)) 0)
(eig-kpoint k-point) (eig-kpoint k-point)
(eig-match-freq? true) (eig-match-freq? true)
Line 78: Line 78:
(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) s)) (* 1.2 sw) 0))))+ (size (* 1.2 (+ (* 2 a) s)) (* 1.2 a) 0))))
- (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 s) dpad) 0) (size 0 sy))+ (center (* 0.5 s) 0) (size 0 a))
(make force-region (direction X) (weight -1.0) (make force-region (direction X) (weight -1.0)
- (center (+ (* 0.5 s) sw dpad) 0) (size 0 sy))))+ (center (+ (* 0.5 s) a) 0) (size 0 a))))
 + (define-param runtime 5000)
(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 a single flux plane (spanning 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 flux planes will overlap thus complicating the analysis whereas the total power through a single flux plane need, by symmetry, 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 for computing the MST, we chose instead to compute the MST along <i>two</i> <code>y</code>-oriented lines (with different <code>weight</code>s to correctly sum the total force), one on each side of the waveguide since by symmetry we need not consider the force in the <code>y</code> direction. +(<b>NOTE</b> if MPB 1.5 is not installed and/or linked with Meep, simply remove all lines pertaining to <code>changes-sources!</code> and the simulation will still work, though with slightly reduced accuracy.)
-We run this simulation over a range of separation distances and compare the result to that obtained from MPB which shows good agreement. Note the presence of mainly repulsive forces for the anti-symmetric mode and mainly attractive forces for the symmetric mode.+There are two important items to note here. The first is that we have defined a <i>single</i> flux surface (spanning 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 flux planes will overlap thus slightly complicating the analysis whereas the total power through a single flux plane need, by symmetry, only be halved in order to determine the value for just one of the two waveguides. The second point to note is that instead of defining a closed, four-sided "box" surrounding the waveguide for computing the MST, we chose instead to compute the MST along just the two <code>y</code>-oriented sides with different <code>weight</code>s to correctly sum the total force, one on each side of the waveguide since by symmetry we need not consider the force in the <code>y</code> direction. Also, we must mention that the MST computed by Meep does not include an often-used factor of 1/2 which needs to be added in post processing in order to properly compare the results to that from the expression above. Choosing a suitable <code>runtime</code> requires some care since a large value is necessary to obtain a steady-state response but this will also lead to large Fourier-transformed fields, used to compute both the flux and the MST, which might cause instabilities in floating-point arithmetic.
 + 
 +We run this simulation over a range of non-zero separation distances and compare the result to that obtained from MPB which shows good agreement. The MPB data for this plot was generated using this [http://ab-initio.mit.edu/~oskooi/wiki_data/parallel-wvgs-mpb.ctl control file] and [http://ab-initio.mit.edu/~oskooi/wiki_data/run_wvgs_mpb.sh shell script]. The plot of the MPB data was generated using this [http://ab-initio.mit.edu/~oskooi/wiki_data/MPB_data_plot.ipynb iPython notebook].
[[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.]]

Current revision

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 our mode-solver package MPB. Note that this will require Meep 1.2+ and MPB 1.5+. Our example consists of two identical dielectric waveguides made of lossless Silicon with a square cross section supporting identical modes and separated in air. Due to the symmetry and orientation of the waveguides, 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 couple more and more 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 for this system is that the sign of the force can be tuned to be either attractive or repulsive depending on the relative phase of the two modes. We will demonstrate this effect 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. For more details, refer to Povinelli et al.. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. This expression has been shown to be mathematically equivalent to the MST (see Rakich et al.) 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 the tricky units altogether. Since the total power transmitted through the waveguide is P = vgU / L where vg is the group velocity, L is the waveguide length and U is defined as before, we focus instead on the force per unit length and power (F / L)(ac / P) where a is an arbitrary unit length and c is the speed of light. This dimensionless quantity enables us to compute both the flux and the force in a single simulation.

We can therefore compute the optical gradient force in two ways, one indirect and the other direct: 1) using any mode solver (in our case MPB), we compute the eigenmode frequency (and corresponding 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 both the force using the MST and 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's set up the two-dimensional computational cell (although since we are interested in a propagating mode of the structure which requires an axial wavevector, this will technically be a three-dimensional simulation)

 (set-param! resolution 30)
 (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 a 1.0) ; waveguide width
 (define-param s 1.0)  ; waveguide separation distance
 (set! geometry (list 
       (make block (center (* -0.5 (+ s a)) 0)
             (size a a infinity) (material Si))
       (make block (center (*  0.5 (+ s a)) 0)
             (size a a 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)
 (set! symmetries (list 
        (make mirror-sym (direction X) (phase (if xodd? -1 +1)))
        (make mirror-sym (direction Y) (phase -1))))

Next, we set the Bloch-periodic boundary condition in order to excite a 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 apriori what the eigenmode frequency will be at any given separation distance, we first excite a spectrum of frequencies using a broadband volume source in each waveguide and then determine the resonant frequency using harminv. The use of a Bloch periodic boundary condition in the z direction will mean that the excited mode (if any) will propagate indefinitely in time which is why we stop the simulation at a fixed 200 time units after the Gaussian-pulsed sources have turned off.

 (define-param fcen 0.22)
 (define-param df 0.06)
 (set! sources (list 
         (make source (src (make gaussian-src (frequency fcen) (fwidth df))) 
                (component Ey) (center (* -0.5 (+ s a)) 0) (size a a 0))
         (make source (src (make gaussian-src (frequency fcen) (fwidth df))) 
                (component Ey) (center (* 0.5 (+ s a)) 0) (size a a 0) 
                (amplitude (if xodd? -1.0 1.0)))))
 (run-sources+ 200 
     (after-sources (harminv Ey (vector3 (* 0.5 (+ s a)) 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 then use this value to accurately excite the mode of interest using the eigenmode-source feature and also compute the force on the waveguide and the power in the guided mode only at this value. The eigenmode-mode feature invokes MPB via a library routine in order to compute the relevant mode of interest and subsequently imports the (steady-state) field profile into the Meep simulation for use as the initial amplitude of the source. This enables an efficient excitation of the relevant eigenmode of interest in the time domain to a much higher degree of accuracy than would otherwise be possible had we simply used a point-dipole source. For more details, refer to chapter 4 of our book.

 (reset-meep)
 (change-sources! (list
    (make eigenmode-source
      (src (make gaussian-src (frequency f) (fwidth df)))
      (component Ey)
      (size a a 0)
      (center (* -0.5 (+ s a)) 0)
      (eig-kpoint k-point)
      (eig-match-freq? true)
      (eig-parity ODD-Y))
    (make eigenmode-source
      (src (make gaussian-src (frequency f) (fwidth df)))
      (component Ey)
      (size a a 0)
      (center (* 0.5 (+ s a)) 0)
      (eig-kpoint k-point)
      (eig-match-freq? true)
      (eig-parity ODD-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 a) s)) (* 1.2 a) 0))))
 (define wvg-force (add-force f 0 1
   (make force-region (direction X) (weight +1.0) 
         (center (* 0.5 s) 0) (size 0 a))
   (make force-region (direction X) (weight -1.0) 
         (center (+ (* 0.5 s) a) 0) (size 0 a))))
 (define-param runtime 5000)
 (run-sources+ runtime)
 (display-fluxes wvg-pwr)
 (display-forces wvg-force)

(NOTE if MPB 1.5 is not installed and/or linked with Meep, simply remove all lines pertaining to changes-sources! and the simulation will still work, though with slightly reduced accuracy.)

There are two important items to note here. The first is that we have defined a single flux surface (spanning 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 flux planes will overlap thus slightly complicating the analysis whereas the total power through a single flux plane need, by symmetry, only be halved in order to determine the value for just one of the two waveguides. The second point to note is that instead of defining a closed, four-sided "box" surrounding the waveguide for computing the MST, we chose instead to compute the MST along just the two y-oriented sides with different weights to correctly sum the total force, one on each side of the waveguide since by symmetry we need not consider the force in the y direction. Also, we must mention that the MST computed by Meep does not include an often-used factor of 1/2 which needs to be added in post processing in order to properly compare the results to that from the expression above. Choosing a suitable runtime requires some care since a large value is necessary to obtain a steady-state response but this will also lead to large Fourier-transformed fields, used to compute both the flux and the MST, which might cause instabilities in floating-point arithmetic.

We run this simulation over a range of non-zero separation distances and compare the result to that obtained from MPB which shows good agreement. The MPB data for this plot was generated using this control file and shell script. The plot of the MPB data was generated using this iPython notebook.

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