E field amplitudes of PhC/Metal devices

#1

I am simulating the reflection and transmission properties of a 1D photonic crystal with a gold film on top. I am also interested in the decaying electric fields across the device. I have attached my file here. toplighttest_lumerical.fsp (1.2 MB) I have a couple of questions listed below.

1. In my 1D field monitor, the fields become slanted and the amplitude is not 1 (even though the plane wave source amplitude is set to 1). I believe this is due to the monitor picking up the plane wave source and reflected light, but I want to double check to make sure it is not an error. I took the device out and ran the sim, resulting in an E field amplitude of 1 ±1%. Shouldn’t this be exact? I noticed that the mesh override region made a small difference in the results.
2. When I switch the background index from air to water (RI=1.335), the resonance peak only changes in intensity, not in spectral location (opposite of what I expected). Could this indicate a faulty setup somewhere?
Thanks for the help!

#2

1. You are correct, the monitor see both the incoming light and the reflection, and the pattern you see is most likely an interference pattern. The variations you observed (±1%) can just be numerical noise. There’s different sources of error in the FDTD algorithm, that can affect the calculation. For instance, I noticed some ripples in the spectra, that seem to be caused by reflections from the PML. Switching to the steep angle profile seemed to get rid of these ripples. You can find some information about these sources of error in the convergence testing page.
2. I’m not too sure about this one. It would be interesting to compare the results of the FDTD calculation with an analytical calculation. For instance, you can try and use the stackrt command to calculate the reflection and transmission spectra.

There might be a slight shift, but it would be very small indeed.

#3

Thank you of the update. How successful were you at getting rid of the ripples and how many PML layers did it take? I ask because more rippling occurs as the Au film thickness increases.

#4

To get the plot I used the default steep angle profile, with 12 layers. That said, you can increase the number of layers to see if it improves the results.

To set these parameter with the script, you can use:

set("pml profile", 3);
set("pml layers", 24);


The pml profile can be 1 (standard), 2 (stabilized), 3 (steep angle), 4 (custom).

#5

I tried that, but did not see any significant improvement even for 200 layers. Also, how do I choose the length of the simulation time found in the FDTD domain?
Thanks!

#6

Yeah, I noticed that too, after my post. The simulation time should be long enough for the light to leave the simulation area. The solver calculates the amount of energy left and stops the simulation when it reaches minimum value set in the advanced tab of the FDTD properties:

The default value is 1e-5. Sometimes it can be necessary to reduce that value (for instance, 1e-8).

If the simulation time is too short (the simulation doesn’t reach the auto-shutoff criteria), it can create some ripples as well (due to the fields being too high when doing the Fourier transform on the time domain fields to get the frequency domain results).

In your case, you can try and reduce the auto-shutoff min, and check the simulation runs long enough (you can see that in the log file).

#7

Thanks, I’ll give it a shot. In the mean time, I have been looking at FDTD’s fit for Palik’s Au data compared to the other gold data sets. The # of max coefficients introduces strange dips and peaks between 200 and 1500 THz in the Palik curve fit, but not in the CRC or J&C data (in spite of the data looking approximately the same. Can you explain what is happening and include details on what FDTD Solutions actually uses when the simulation is running?
Thanks!

#8

When using dispersive materials, we have a well known relationship between the displacement field and the electric field, in the frequency domain:

FDTD is a time domain technique, so this relationship becomes a convolution product:

So we need to express the permittivity in a way that:

1. Must be able to solve this convolution product in the time domain
2. Must be causal (it must satisfy the Kramers-Kronig relationships)

This puts some limitations on the forms of the permittivity we can use. Our multi-coefficient material model fits the material data with a set of functions that fulfils these conditions. The fit is what the solver will actually use in the calculation.

Obviously, you want the fit to be as close as possible to the actual data but, as you noticed, it can be difficult to get a good fit. In your case, although the data seem similar, the differences make the fit very poor for the Palik data over the bandwidth of interest. You can try and improve the fit, but keep in mind the data come from experiments and can include noise and errors.