Optical Forces

This tutorial demonstrates Meep's ability to compute classical forces using the Maxwell stress tensor (MST) as well as the eigenmode source feature which involves the mode solver MPB. See parallel-wvgs-force.py. The geometry consists of two identical silicon waveguides with square cross section in vacuum (shown in the figure below). Due to the orientation of the waveguides, the two modes can be chosen to be either symmetric or anti-symmetric with respect to a mirror-symmetry plane between them. As the two waveguides are brought closer and closer together, their modes couple more and more and give rise to an optical gradient force that is transverse to the waveguide axis. This is different from radiation pressure which involves momentum exchange between photons and is longitudinal in nature. An interesting phenomena that occurs for this system is that the force can be tuned to be either attractive or repulsive depending on the relative phase of the modes. We will demonstrate this effect in this tutorial.

The optical gradient force on each waveguide arising from the evanescent coupling of the waveguide modes can also be computed analytically:

where ω is the eigenmode frequency of the coupled-waveguide system, is the separation distance between the parallel waveguides, is the conserved wave vector and is the total energy of the electromagnetic fields. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. For more details, see Optics Letters, Vol. 30, pp. 3042-4, 2005. This expression has been shown to be mathematically equivalent to the MST in Optics Express, Vol. 17, pp. 18116-35, 2009. We will verify this result in this tutorial.

It is convenient to normalize the force so as to eliminate the tricky units altogether. Since the total power transmitted through the waveguide is where is the group velocity, is the waveguide length and is defined as before, we focus instead on the force per unit length and power where is an arbitrary unit length and is the speed of light. This dimensionless quantity enables us to compute both the flux and the force in a single simulation.

We can compute the optical gradient force in two ways: (1) using MPB, we compute the eigenfrequency and group velocity for a given mode at different separation distances and then use a finite-difference scheme to numerically evaluate the formula from above, and (2) using Meep, we compute both the MST and the power transmitted through the waveguide for the guided mode. In this particular example, we consider just the fundamental y-odd mode which shows the bi-directional force.

First, we define the necessary modules:

from __future__ import division

import math
import meep as mp
import argparse

def main(args):

Let's set up the computational cell. The waveguide cross section is 2d but since we are interested in a guided mode of the structure which requires specifying the axial wavevector, this will in fact be a 3d simulation:

resolution = 30
nSi = 3.45
Si = mp.Medium(index=nSi)
dpml = 1.0
sx = 5
sy = 3

cell = mp.Vector3(sx + 2 * dpml, sy + 2 * dpml)
pml_layers = mp.PML(dpml)

a = 1.0     # waveguide width
s = args.s  # waveguide separation distance

geometry = [mp.Block(center=mp.Vector3(-0.5 * (s + a)),
                     size=mp.Vector3(a, a, 1e20),
            mp.Block(center=mp.Vector3(0.5 * (s + a)),
                     size=mp.Vector3(a, a, 1e20),

Two mirror symmetries can be used to reduce the size of the computational cell by a factor of four:

xodd = args.xodd
symmetries = [mp.Mirror(mp.X, phase=-1.0 if xodd else 1.0),
              mp.Mirror(mp.Y, phase=-1.0)]

Next, we set the Bloch-periodic boundary condition for the mode with wavevector π/:

beta = 0.5
k_point = mp.Vector3(z=beta)

Since we do not know apriori what the eigenfrequency is at a given separation distance, we use a broadband source to excite multiple frequencies and then use harminv to determine the resonant frequency. The use of Bloch-periodic boundary conditions in the direction means that the excited mode (if any) will propagate indefinitely in time which is why we stop the simulation at 200 time units after the sources have turned off. This number is somewhat arbitrary.

fcen = 0.22
df = 0.06
sources = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), component=mp.Ey,
                     center=mp.Vector3(-0.5 * (s + a)), size=mp.Vector3(a, a)),
           mp.Source(src=mp.GaussianSource(fcen, fwidth=df), component=mp.Ey,
                     center=mp.Vector3(0.5 * (s + a)), size=mp.Vector3(a, a),
                     amplitude=-1.0 if xodd else 1.0)]

sim = mp.Simulation(resolution=resolution,

h = mp.Harminv(mp.Ey, mp.Vector3(0.5 * (s + a)), fcen, df)

sim.run(mp.after_sources(h), until_after_sources=200)

f = h.modes[0].freq
print("freq:, {}, {}".format(s, f))

Once we have determined the eigenmode frequency, we then excite this mode using the eigenmode-source feature. Also, we compute the force on each waveguide and the power in the guided mode. The eigenmode-mode feature invokes MPB (which needs to be available as a library) to compute the relevant mode of interest. The steady-state field profile from MPB is imported into Meep for use as the initial source amplitude. This enables an efficient excitation of the eigenmode 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 Section 4.2 ("Incident Fields and Equivalent Currents") in Chapter 4 ("Electromagnetic Wave Source Conditions") of the book Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology.


new_sources = [
        mp.EigenModeSource(src=mp.GaussianSource(f, fwidth=df), component=mp.Ey,
                           size=mp.Vector3(a, a), center=mp.Vector3(-0.5 * (s + a)),
                           eig_kpoint=k_point, eig_match_freq=True, eig_parity=mp.ODD_Y),
        mp.EigenModeSource(src=mp.GaussianSource(f, fwidth=df), component=mp.Ey,
                           size=mp.Vector3(a, a), center=mp.Vector3(0.5 * (s + a)),
                           eig_kpoint=k_point, eig_match_freq=True, eig_parity=mp.ODD_Y,
                           amplitude=-1.0 if xodd else 1.0)


flx_reg = mp.FluxRegion(direction=mp.Z, center=mp.Vector3(), size=mp.Vector3(1.2 * (2 * a + s), 1.2 * a))
wvg_pwr = sim.add_flux(f, 0, 1, flx_reg)

frc_reg1 = mp.ForceRegion(mp.Vector3(0.5 * s), direction=mp.X, weight=1.0, size=mp.Vector3(y=a))
frc_reg2 = mp.ForceRegion(mp.Vector3(0.5 * s + a), direction=mp.X, weight=-1.0, size=mp.Vector3(y=a))
wvg_force = sim.add_force(f, 0, 1, frc_reg1, frc_reg2)

runtime = 5000

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('-xodd', type=bool, default=True, help='odd mirror-symmetry plane in the X direction?')
    parser.add_argument('-s', type=float, default=1.0, help='waveguide separation distance')
    args = parser.parse_args()

There are two important items to note: (1) We have defined a single flux surface to compute the Poynting vector along which spans an area slightly larger than both waveguides rather than two separate flux surfaces, one for each waveguide. This is because in the limit of small separation, two flux surfaces overlap whereas the total power through a single flux surface need, by symmetry, only be halved in order to determine the value for just one of the two waveguides. (2) Instead of defining a closed, four-sided "box" surrounding the waveguides for computing the MST, we chose instead to compute the MST along two -oriented sides with different weight values to correctly sum the total force one on each side of the waveguide. By symmetry, we need not consider the force along the direction. Choosing a suitable runtime requires some care. A large runtime is necessary to obtain a steady-state response but this will also lead to large values for the discrete Fourier-transformed fields used to compute both the flux and the MST. These large values might may produce floating-point errors.

We run this simulation over a range of non-zero separation distances and compare the result to that obtained from MPB. This is shown in the figure below. The two methods show good agreement. The MPB data for this plot was generated using this Scheme file and shell script. The plot of the MPB data was generated using this Jupyter notebook (html).