How to calculate Optical Forces using MEEP

The optical force acting on an arbitrary particle in an arbitrary incident field can be calculated using the integral

$\mathbf{F} = \int_{S} \left\lbrace \mathbf{T}(\mathbf{r},t) \right\rbrace \cdot \mathbf{n}(\mathbf{r}) dS$

where $\left\lbrace \cdot \right\rbrace$, $S$ and $\mathbf{T}$ represent the time-averaged operator, surface enclosing the particle and Maxwell's stress tensor, respectively. Maxwell's stress tensor is given as

$\mathbf{T} = \left[ \varepsilon_0 \mathbf{E}\mathbf{E} - \mu_0 \mathbf{H}\mathbf{H} - \frac{1}{2}(\varepsilon_0 E^{2} + \mu_0 H^2) \mathbf{I} \right]$

which requires all field components to be calculated. You can obtain the solved field components in a number of ways. We can use the data from MEEP to calculate the stress tensor and hence the optical forces. If we are interested in finding the forces at a grid of points then we must run the MEEP simulation for each point.

I have provided a set of files which calculates the forces a particle experiences in a standing wave. It should be fairly easy to modify to suit your simulation.

- calculateforce.py : this file runs everything and plots the forces. You will need numpy/scipy installed.

- standing.ctl : meep simulation file that generates stress tensor for a particle at different locations (you can see txx, txy, etc definitions)

- standing-full.ctl : meep simulation to obtain field density (only for plotting)

- h5tonumpy : convert h5 dataset to numpy file.

This method of obtaining optical forces is very slow and there are inherent inaccuracies. There was a paper that mentioned this. I will provide I link to it when I find it again. Essentially we are only using FDTD to obtain the fields which can be obtained using other scattering codes/methods. However, for more complicated geometries and/or materials (such as plasmonics/meta-materials) the FDTD brute force method can be useful.

Please note: The latest version of MEEP now has an inbuilt function to calculate Maxwell's stress tensor. It is much easier to calculate forces using that. However I leave this here for legacy reasons.