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

where , and represent the time-averaged operator, surface enclosing the particle and Maxwell's stress tensor, respectively. Maxwell's stress tensor is given as

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.