Problem with Kerr NL materials in simulation


im trying to simulate a normal plane wave in a kerr medium by using a Chi3 material.
Even though my source is a plane wave, it gets distorted as it propagates, i think this is a problem with the Kerr material model because when i use a linear material instead it does not happen.

I would greatly appreciate it if you could help me with this problem.
temp2.fsp (287.5 KB)

Reflection from PML boundary condition
made this topic public #2


Hi @nimrodg

Thank you for reaching out.

Can you please elaborate more on what you mean by distortion on the plane wave source? Here is a plot of propagated mode with nonlinear material extracted from time monitors, but I am not quite sure of the problem (maybe something very obvious that I am missing?):

Sorry that I was not much of help, but please let me know of your thoughts.


Yes, the propagating wave in the simulation looks like that, and that should not happen.
A plane wave propagating in a kerr medium should remain a plane wave with constant field amplitude.
In the simulation you can see that the amplitude of the plane wave has growing fluctuations in it.

Im not sure why this happens, my only guess is that maybe the kerr model in Lumerical changes the refractive index according to real(E) and not abs(E). I tried looking at the code in the Chi3 medium but it did not clarify anything.
Any thoughts?



Dear @nimrodg

Thank you for providing additional information. I understand the problem now.

I spend quite a lot of time with a colleague of mine on this case. We tried a couple of different scenarios, and we uncover some strange behaviour that require a deeper look in the model with our R&D team. But in the mean time, I guess I found a solution/work around.

First, I added a point monitor to the center of simulation region. You can use this point monitor to look at the electric field spectrum. While Kerr nonlinearity is meant to affect only the optical refractive index, for fine x-mesh size of 5 nm, I observed some higher frequency components:

temp2_modified.fsp (288.8 KB)

These unwanted frequencies cause beating in time domain signal and do not exist for linear material. I repeated simulations for different types of plane wave source using the CW source script, but could not obtain any better/convincing results.

Finally, I noticed that the higher frequency components can be removed using a coarser mesh. For example, when I used a mesh size of 50nm, I obtained a signal with almost flat amplitude:

temp2_modified_ver2.fsp (1.2 MB)

Note that this is slightly different simulation file than previous one and I used the script below to create a time signal:

# this script shows how to create a custom time signal

# specify signal

  t        = linspace(0,150000e-15, 40000);
  lc = 500e-9; # wavelength
  f = c/lc;
  w_center = 2*pi*f; 
  delta_t  = 100e-15;
  offset1   = 500e-15;

  # define amplitude and phase
  amp = (t < offset1)*exp(-0.5*(t-offset1)^2/delta_t^2) + (t >=offset1);
  phase = -w_center*t;

  # set the signal for the desired source

As I mentioned, I will need to dig more into the problem to see what was causing those additional frequencies. However, it looks like using a time monitor to look at the spectrum and a coarser mesh to suppress high frequency oscillations will be a good plan.

I hope this was helpful, and sorry for any inconvenience.

Please keep me updated if you have any additional questions.



OK, thats good to know.
Im afraid the coarser mesh might hurt other results of the simulation but ill try it.
Is there any specific reason why you wrote your own plane wave script rather than use the default one? Should i do that too?
I also wanted to ask how exactly does your Kerr material work, does it use the RMS value of the electric field or the absolute value on order to calculate the change of refractive index?

Thank you for the help, and please let me know if there are any updates about this issue.


I have another request if possible.
Can you put me in touch with one of the designers\programmers who is dealing with this issue so that i can talk to them directly? It would help me a lot to understand the mechanism that is used to get the Kerr effect.

Thank you.


Dear @nimrodg

After discussing your case further with my colleagues, I think we finally found the source of the problem.

In your simulation file, from the advanced options tab of Edit FDTD windows you had “always use complex fields” box checked. This option should be unchecked for all nonlinear simulations:

Quoted from our internal discussions:

It will definitely create problems when using complex fields and nonlinear simulation. Complex fields in the time domain are not physical, but in linear systems it is irrelevant. In nonlinear systems it is a bit unclear what to do: does E(t)^2 mean |E(t)|^2 or something else?

I was also wrong to say that Kerr nonlinearity is meant to affect only the optical refractive index. In fact, you expect to create all the harmonics at (2m+1)*w-source with Kerr material. However, using complex fields suppresses half of the harmonics, while working with the real fields will return all the harmonics.

Quoted from our internal discussions:

The plugins will be called twice with complex fields – once for the real part of E and again for the imaginary part of E. It is interesting that this suppresses half the harmonics – I will have to think that through at some point.

If always use complex fields box is unchecked, your simulations will include only the real part of time domain electric field (i.e. E*cos(wt)), and thus you expect the time field to oscillate in time. If you are interested to obtain the envelope, you can use the script below which basically filters higher order frequencies beyond Nyquist limit:

mname = "monitor";

t = getdata(mname,"t");
Ey = getdata(mname,"Ey");
Ey = pinch(Ey);

# fourier transform
Eyf = fft(Ey);
w = fftw(t);

# apply filters in frequency domain
w1 = 2*pi*c/500e-9;
#w2 = 2*pi*c/1500e-9;

filter_complex = 2*((1:length(w)) <= length(w)/2);
filter1 = 2*exp(-(w-w1)^2/(150e12)^2);
#filter2 = 2*exp(-(w-w2)^2/(150e12)^2);

# plot the spectrum and filters
plot(w/(2*pi)*1e-12,abs(Eyf)/max(abs(Eyf)),filter1,filter_complex,"f (THz)","spectrum and filters");
legend("spectrum","filter1","filter2","filter complex");

# calculate the pump and lasing field in the time domain
Ey_complex = invfft(Eyf*filter_complex);
Ey_w0 = invfft(Eyf*filter1);
#Ey_2w0 = invfft(Eyf*filter2);
# get new time vector, different due to zero padding
t2 = fftw(w);

# plot the pump and lasing output
plot(t2*1e12,abs(Ey_w0),"t (ps)", "|E|");

plotxy(t*1e12,abs(Ey),t2*1e12,abs(Ey_complex),"t (ps)", "|E|");
legend("|E(t) original|","envelope");

While this should fix the problem for you, below I have tried to answer your questions:

I used a time domain signal with a slow increase in amplitude. For pulses that are truncated or have a sudden start, you might create some other frequencies when you Fourier transform them. This is something to bare in mind if you saw some other non-sense frequencies.

You can use finer mesh now. If you go to a really small mesh, you can get more reflections from the pml, maybe increase the number of layers to 24 or 32 so you don’t have to worry about it. I don’t think this is the problem here, but it is a good idea.

We invite our developers to respond directly to some posts, but I think it will be faster if I stay in charge at the moment as our developers might be quite busy before big releases. But please let me know if you had further questions and I will discuss it further with our team.

I hope this answered your questions.

Nonlinear material plugin with Bloch boundary conditions - four wave mixing application example
Plane wave source size changes
Reflection from PML boundary condition
Isotropic nonlinear material