Optical force

Hi There, My first experience with Lumerical! I appreciate your help in these…

1)when we define a material ourselves, with n and k, to show that the
material is lossy shall I take k a positive number or negative number?

  1. how can I calculate the field, field gradient, and probably any other
    function of the field values calculate by Lumerical at certain point in
    space within the software.

3)I am attaching this file where I am trying to calculate the optical
force on the particle and also need the field and field gradient
components at the center of the particle.
… Well just don’t see where I can attach my files.But hopefully you can help me with that.


To model a lossy material, the value of the extinction coefficient, k should be positive.

You can use different monitors in FDTD and MODE solutions to capture the field profile at any point. Once you have the field profile, you can use the script environment to calculate secondary quantities such as the gradient of the field. To learn more about the different monitors take a look at this KB page.

You can click the “upload” button to upload your file when you post in the forum (see screenshot below).


For calculating the optical force (or the Maxwell Stress tensor), Lumerical has a knowledge base article.


With a simulation file plus script, is a place to start. Also, you can add Maxwell Stress tensor (either on the box or volumetric) from the advanced monitor tab.



I would like to use the Maxwell Stress tensor’s script to calculate the optical force on the waveguide’s surface due to the guided modes. Could you give any tips ?




If I understand correctly, you would like to simulate something similar to the following publication:

If so, you should be able to get the net force on the waveguide by placing the optical force analysis group so that the monitors enclose the cross section of the waveguide but does not enclose the substrate below or any other structures surrounding the waveguide. Then by dividing the calculated force by the span of the analysis group along the waveguide propagation direction you can get the force per unit length on the waveguide.


  You are rigth, it is this kind of force that I am talking about. More precisely, I am trying to reproduce the results from Fig. 3 presented in: [https://www.researchgate.net/publication/49782199_Scaling_of_optical_forces_in_dielectric_waveguides_Rigorous_connection_between_radiation_pressure_and_dispersion] or  [https://www.osapublishing.org/ol/abstract.cfm?uri=ol-36-2-217] . 

  Since my first post I have tried to do exactly what you said. However, even dividing by the waveguide length I am off by several orders of magnitude. Please, could you take a quick look in my fdtd file ?        




Unfortunately I did not have access to the publication that you referred to, but if you attach a copy of your simulation file, I could definitely take a look.

Here is the link:

Please, let me know as soon as you download these files, so I can remove the link.

Thank you very much.


I’ve just downloaded the files, so you can go ahead and remove the link it you’d like.

After a quick check, I have a couple of suggestions to make:

  1. I would recommend extending the x max position of the waveguide into the x max boundary of the simulation region. This way, the PML boundary will absorb the incident light. With the current setup, once the injected mode reaches the end of the waveguide, some of the light will be reflected, forming a cavity where the light bounces back and forth inside the waveguide.
  2. From the diagram on figure 3 of the publication, it looks like px and py refer to the pressure on one surface of the waveguide. In this case, you can set the span and position of the optical force analysis group to enclose just one the surfaces at a time. For example, to get py, you can have the group enclose only half of the waveguide from y=0 to y=0.2 um and extract the force in the y-direction. Otherwise, if the mode is symmetric and you include the full waveguide structure, the net force on the object which includes both the y min and y max interface of the waveguide will cancel out.

Please try it out and let me know if you’re able to get similar results after these changes!

1 Like

Thank you very much. I extended the waveguide length to the PML region and I also set the span and the position of the optical force analysis group to enclose only the right hand side of y-oriented waveguide surface in order to get only py (which is the main one, as you can see from the force density distribution showed in Fig. 1 (f) of the publication). By doing these changes I have gotten a force of 1,25E-24 N, and dividing it by the span of the analysis group along the waveguide propagation direction (1 um in x-direction) I got a py of 1,25E-18 N/m, which still out by a factor of E10 from the publication.

I looked into the simulation setup in more detail, and one thing I noticed is that the mode which was being injected in the simulation file that you provided corresponded to the TM mode instead of the TE mode, and the y and z span of the simulation region could be extended a bit to make sure that the mode profile isn’t truncated at the boundaries of the simulation region. However, I tried making these changes and the result still isn’t very close to the reported result in the publication.

At this point, I’m not sure what the main problem is, but I may be able to look into the calculation in more detail and get back to you early next week.

Thank you very much for your help so far. I’m also trying to figure out what is happening.

I’ve taken another look and I had made a mistake in the normalization of the results last time. With the following simulation file I’m able to get a result which is on the order of the force which is reported in the paper:
modified_MST_in_a_waveguide.fsp (1.2 MB)

I used the following script to extract the normalized results so that the results are in the units reported by the paper:

# frequency
f = c/1.55e-6;

# power injected by source in Watts
sp = sourcepower(f);

# area of y-normal surface where optical force object is placed over in m^2
A = 0.315e-6*0.5e-6;

# get y-component of force from analysis group in N
F_total_dataset = getresult("optical_force_MST","F_total");
F_total = F_total_dataset.F_total;
Fy = F_total(2); # y component is the second component of the matrix

# normalize result to units of N/m^2/W and print to script prompt
F_norm = Fy/A/sp;
?"Normalized force in y direction is: "+num2str(F_norm);

Please try it out!

1 Like

Thank you very much. I will try it.