PML reflection that is prevalent on a specific wavelength range

I am using Lumerical FDTD Solutions 8.15.786 on Windows 7 64-bit. I have a geometry that contains a semi-infinite slab of air on top of a semi-infinite slab of silver with a slit perforated on it. There is air inside the slit and I use a TFSF source that spans a wavelength range from 400 nm to 1200 nm. All the boundaries are PML and I pay particular attention to have some space between TFSF boundaries and the simulation domain boundaries. On top of that I use a uniform simulation mesh in all directions without any mesh-override region. The TF/SF is excited from the air side, and it propagates down the slit perforated on the infinite metal slab.

My question is about how to implement the correct PML settings so that there are no reflections from the bottom of the metal, since I try to simulate a case with a metal of infinite thickness. However, when I use the modal fields to calculate the forward and backward propagating mode amplitudes, I found out that there is also substantial reflection from the bottom of the metal in a wavelength range between 850 nm to 1200 nm.

I tried out various approaches to remedy this issue: increasing the number of PML layers, changing the PML type from standardized to stabilized since there are material boundaries on the PML. However, these solutions did not work and in one of the cases(I think it was the one using standard PML) I got a diverging simulation warning even though I increased the number of PML layers. I use the streched coordinate PML that is available in the new releases of this software.

The simulation geometry in the YZ cross-section and the plot of the unnormalized reflection coefficients against the wavelength are attached as images. In these attached graphs A’s are the amplitudes of the forwards propagating modes, whereas B’s are amplitudes of backward propagating modes. I would appreciate it if you could comment on why there is reflection and how to alleviate it.

Plot of A(forward) B(backward) propagating mode coefficients

EDIT: The simulation files are attached to this post.
Slit_f2a_PML7.fsp (256.4 KB) This file is the latest simulation configuration I use
Slit_f2a_PML1.fsp (250.2 KB) This file is the original file I sent to support


I have some comments and suggestions (see modified file: metal_slit - 0p9-1.fsp (250.0 KB)):

  1. You can use the TFSF source but you need to be aware that the transmission measured by power monitors will be normalized to the source power, which depends on the injection area of the TFSF source. You could compare the transmission at different places along the slit or compare the results between different slit sizes, and the results will be meaningful because the comparison (for example the ratio between two transmission values) won’t be affected by the normalization issue. Alternatively, if you are trying to simulate a beam that is incident on the slit you could use a Gaussian source, which also has a finite width. At the end of the day, I would say the shape of the fields in the slit won’t depend much on the way you excite them.
  2. The issue with the back reflections of the PML usually happens when you have a dispersive material that extends through the PML. The stabilized profile for the PML where this happens is the first thing you can try, but as you discovered it can lead to artifacts in the transmission. These artifacts can be minimized by increasing the number of layers; however, I found that even if you do that the auto shutoff level does not decrease enough (it goes down to 1e-3, stays there for a long time and then slowly increases), so the frequency-domain results might not be correct.
  3. I tried running a simulation for a smaller wavelength range (from 0.9 to 1um) instead of the full range you were using (0.4 to 1.2um). In this range the simulation ran nicely and the autoshutoff was triggered. To check if there were any back reflections from the PML I placed a mode expansion monitor inside the slit. This monitor calculates how much forward and backward transmission you get into the mode supported by the slit (I just looked at the fundamental mode, which is the one that I expect to be excited more strongly):
    Note that T_forward matches T_total (the total transmission measured by the power monitor used in the analysis) and there is no backward transmission. I suspect that there is some region in the material fit that leads to absorption problems in the PML. When the bandwidth of the pulse is smaller these problems are not triggered so this might be the reason why using a smaller bandwidth works.
  4. You are probably interested in looking at the whole wavelength range 0.4 to 1.2um. If that is the case, you can run a sweep to cover the whole range with smaller range simulations (like the one I did for 0.9 to 1um). Just make sure you use the “set simulation bandwidth” so that the mesh and material fit is used consistently when sweeping over the complete range:
    This approach is discussed in this posting.
  5. One more reason for using simulations with a smaller wavelength range is that strictly speaking the mode expansion monitor only gives you correct results for the center frequency because the field profile used in the analysis is the one for the center frequency. This is the same issue that happens with mode sources. Therefore, if you want to use the mode expansion monitor to find the forward and backward propagation you should use smaller wavelength ranges. I am not sure what method you used to calculate the back reflection but I think the mode expansion can give you an accurate result.

Hope this helps!


Dear Federico Gomez,

Thanks for your reply. I calculated the modal
fields with Comsol which solves equations with the FEM method in
frequency domain. Also the time convention in Comsol is opposite and I
took the appropriate precaution by conjugation every electromagnetic
field. Regarding the points you raised in your reply, please see my comments below:

1. The fields at the surfaces with plane wave excitation are also needed;
therefore using a TFSF source will be best. For transmissions I
normalize with a transmission monitor in empty space with the same
simulation configuration in my regular calculations and it works fine.
I also tried disabling the extension of dispersive material through the
PML layer and terminating the metal just at the FDTD boundary before PML
layers, yet this also did not remedy the situation. Note that this approach was suggested in Lumerical knowledge base.

2. I also tried the mode expansion monitors and inspect the results with my own
calculation method to compare the two. Essentially I have a similar
mathematical approach with the mode expansion monitors and I think at
least the trend in magnitude of A and B should be similar in both calculations. In one of my finite simulations the trend was similar in some cases, yet the magnitudes were different due to normalization.

**3.**Is it easy to merge all the results obtained with different simulations over small wavelength
ranges? Can the merge be done inside FDTD Solutions?

**4.**This also seems to explain why I got a discrepancy between my calculations
and the A, Bs obtained via the mode expansion monitors in some cases.

Furthermore, I noticed that upon shrinking the simulation domain in X and Y I get different results in the mode expansion monitors compared to the simulation having a larger domain in X and Y. I would really appreciate suggestions on how to eliminate back reflection and obtain consistent results over the broad wavelength range I am interested in.

P.S: Dear Federico the original file I sent for inspection was Slit_f2a_PML1.fsp the other file, namely Slit_f2a_PML7.fsp contains a shrunk domain in X and Y, yet the range in Z is increased. This seems to change the electromagnetic fields in cuts with a Z-normal drastically hinting that PMLs on the sides might also be problematic.

Do the settings eliminate discontinuities and optimize for short pulse can have anything to do with these problems? By the way how can I run a parameter sweep and save all the files corresponding to different wavelength ranges for access later on?

Dear Ongun,

I had not tried disabling the “optimize for short pulse” option. When you disable it, the simulation converges for the complete wavelength range 0.4 to 1.2um! This is another indication that the problem that was causing the divergence was some feature in the material fit outside the desired range because the power injected for wavelengths outside the desired range is minimized when the “optimize for short pulse” option is disabled. These are the results for the full wavelength range using the mode expansion monitor in metal_slit - 0p9-1_TFSFsource.fsp (250.2 KB):

It seems like there is some back reflection at wavelengths close to 0.9um. However, it is important to remember that the mode expansion monitor only uses the mode calculated at the center frequency so the results for other frequencies might not be as accurate. Therefore, I also created a sweep over “single-frequency” simulations (frequency span of the source set to zero) as shown in the attached file: metal_slit - 0p9-1_TFSFsource - sweep.fsp (251.7 KB). The results for the same wavelength points do not show any back reflection and the total transmission agrees very well with the broadband simulation:

Therefore, the back reflection seen in the mode expansion monitor of the broadband simulation is probably just an artifact of the mode expansion and you can trust the total transmission calculated in the broadband simulation.

For tips on sweeps over single-frequency simulations take a look at this thread:

Note that I switched back to the TFSF source, which is fine as long as you are aware of the normalization issues.

Regarding the x and y span of the simulation region it is important to make sure there is a distance of at least half the maximum wavelength between the structure (the slit in this case) and the PML so that the evanescent fields at the PML are minimized. This is probably the reason why you were seeing some differences when you made the x and y span smaller. In the simulation files I attached, I made sure the span is large enough.

1 Like

By the way I use a uniform mesh since it is easier to deal with the dot
products I have in my analysis code. Will this make a difference? My dx is
5 nm.

In my analysis I use normalized fields from the Comsol FEM software RF
module and these fields are in frequency domain. I will double check the
signs of these as Lumerical uses negative time exponential, whereas Comsol
uses positive time exponential for forward propagating waves. I have a few
more questions:

1. How far away from the PMLs should the TFSF boundaries be? Can I have the touching each other without invading the grayed region around the TFSF source? Please have a look at the attached simulation files. Regarding the previous discussion the setup PML34 is correct as I leave enough space around the slit. Is this reasoning correct?

Slit_f2a_PML36.fsp (266.2 KB)
Slit_f2a_PML34.fsp (269.0 KB)

2. In my opinion the simulation time also has an effect on the outcome, so will I benefit
from increasing the simulation time? Is using the Auto Shuttoff feature always justified or is it just safer to wait for the complete simulation to finish for more trustworthy results?

3. Does the slit perforating through the metal a big problem from standard PML
profile? I ask this since the epsilon changes considerably and in your fsp files you use steep angle PMLs apart from the stabilized PML at the z_max. What is the justification for that? Furthermore, can checking extent structure through pml and extending it also manually(By making the metal, slit sections larger and longer, respectively) cause problems?

4. I also plotted the square of the electric field norm along a line in through the center of the slit and noticed that the electric field norm is much higher at the metal PML which is really far away from the injection
plane. I clicked record data in pml and extended the monitor through PMLs to make sure I caught everything. I think if I do not click this option the data after z=600 nm(the boundary of the FDTD region) will not be available. Does this setting work like this? Here is a sample plot:

Hi Ongun,

Sorry for the wait. See some suggestions and comments below:

It should be OK to use the uniform mesh. Do you need it everywhere in the simulation region? There are two main advantages of using the nonuniform mesh:

  1. You save memory and reduce execution time.
  2. You can make the mesh near the PML coarser, which increases the thickness of the PML layers so you get better absorption.

Regarding our phase convention, FDTD solutions uses exp(-iwt+ikz) for forward propagating plane waves.

The position of the TFSF boundaries in your simulation is fine. The important thing is to make sure the side boundaries are not placed beyond the simulation boundaries to avoid the edge effects that you get when using a plane wave with finite size.

Regarding the distance from the slit to the PML it should be at least half the maximum wavelength. In the PML34 file you are leaving more than that distance, which is fine but you might want to reduce it to speed up the simulation.

The autoshutoff feature helps you save execution time, since there is no point in keeping a simulation running when the fields in the simulation region have left or been completely absorbed.

The right time to end the simulation depends on the situation. For most cases, using the default autoshutoff level min (1e-5) works well. However, since this is a kind of average of the energy left in the simulation region, it might end the simulation before some energy localized somewhere in the structure has decayed enough. A typical signal of an early termination of the simulation is finding ripples in the frequency-domain results.

You can always check if increasing the simulation time and reducing the autoshutoff min changes your results. You can then identify the appropriate autoshutoff level that does not cause artifacts in the results while saving execution time.

I just noticed that you actually don’t need the “steep angle” and “stabilized” PML profiles. You can simply use that “standard” one. I was trying the other profiles to avoid the divergence that you were getting initially, but in this case disabling the “optimize for short pulse” option in the source settings is enough to avoid the divergence.

You don’t have to extend the structure manually through the PML if you check the “extend structure through PML” option and make sure the slit and the metal are touching the PML boundary. Just to be safe I always extend the structure through the PML manually.

You are right, that is the way the “record data in pml” works. If it is unchecked (default) the monitor data is cut off by the simulation region boundary.

Are you plotting the field intensity in log scale?

I would say it is normal to have a high field intensity near the upper PML because the fields are guided by the slit into the PML. You can see that they decrease as they get into the PML, as expected.

1 Like

Hello Federico Gomez,

I think I can use uniform mesh only inside the regions where I have the monitors placed. I do not need them in the PML regions as I will not post-process data from there. I think I can make part of the mesh uniform and part of it nonuniform.

Since the phase convention is different(minus that of Lumerical’s), I got the conjugate of each field resulting from COMSOL simulations. I think this will be sufficient to convert from one phase convention to another.

As far as I understand it is as long as the greyed boundaries around the TFSF source does not touch the PML boundaries we are fine. I think the configuration in the simulation PML34 depicts the limit that these boundaries can be pushed the most.

I will go without auto shutoff as I find it easier and I will change all the PMLs to “standard” PMLs. By the way does the last boundary having the slit metal interface not require the stabilized PML?

Okay I just extend it anyways just to make sure that the PML is the same type of material as of the simulation domain.

Indeed, in that plot I actually plotted the electric field norm squared(intensity) in log scale. I also noticed that they decrease, but the PML is very dispersive and this decrease is wavelength dependent. This plot also indicates that there is some reflection around 900 nm.

I think for saving simulation time considerably I can make use of symmetries PEC, PMC boundaries as I am only interested in the first mode for the time being.

This problem is still persistent and it seems to be related with problems in PML absorption around the cutoff wavelength of the waveguide. This is mentioned in the article An Auxiliary Differential Equation Formulation for the Complex-Frequency Shifted PMLby Stephen D. Gedney and Bo Zhao. This article is actually the second reference on the PML boundary conditions page on Lumerical Knowledgebase. However, there is some normalization utilized details of which is omitted as can be seen on the following screenshot.

The same article also mentions that using parameters as follows sigma_max = sigma_optimum, kappa_max = 15(both with polynomial scaling m =3) and a_max = 0.4 which is linearly scaled. . The sigma_optimum is defined in the book Computational Electrodynamics: The Finite-Difference Time-Domain Method by Allen Taflove as follows:

The problem here is with the correct values for a_max and sigma_max since they are normalized according to Lumerical’s convention. If you could provide information on this I think I will be able to remedy the problem

Dear Ongun,

These are some comments and suggestions:

Yes, you can use standard PML for the upper z boundary. However, I would increase the number of layers as explained below.

I see your point. In this case I would suggest increasing the number of PML layers to make sure the fields have decayed enough by the time they reach the outer limit of the PML. I increased the number of layers from 8 to 32 for the zmax PML boundary in this simulation file: metal_slit - 0p9-1_TFSFsource_v2.fsp (254.5 KB). The field intensity is still strong for wavelengths in the range ~0.9 to 1 um, but it decays as you get further into the PML. Also, the mode expansion monitor now shows no significant back reflection.

You can definitely use symmetric/antisymmetric boundary conditions to speed up the simulation (see attached simulation file). I would not use PEC or PMC boundaries as this would lead to extra reflections from the boundaries.

Regarding the choice of optimal parameters:

  • We use sigma values close to the optimal described in the extract from Taflove’s book.
  • Using Eq. 7.67 for dispersive materials and graded mesh can be tricky because this is a frequency-domain expression so you would have a sigma that depends on frequency. This is difficult to implement in the time domain, so you would have to choose one frequency for calculating sigma and there is no guarantee this would work in a broadband simulation.
  • The choice of PML parameters in the profiles we provide (standard, stabilized, etc.) has been the result of extensive tests to find optimal values that work in most cases. You can find a thorough discussion of how the PML parameters work in this webinar:

Hope this helps!