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.