Graphene-based metamaterials

I am trying to perform similar analysis of SiO2 + graphene stack using FDTD Solutions as in the article Tunable spectral and spatial filters for the mid-infrared based on hyperbolic metamaterials. Mainly I would like to achieve similar results as in Fig. 2.

Fig. 2


I was using the 2D graphene model with 2D FDTD region as shown below. There were 20 layers of combined SiO2 (150 nm) + 2D graphene. To perform analysis I was using the S parameters group from Object Library - with plane wave (Bloch/periodic) as a source. The monitors were placed 5 um from the centre of the structure. The FDTD boundaries were set to PML (16 layers) in x direction and periodic in y direction. Additionally there was mesh override region around the structure with dx = 1 nm.

FDTD Solutions view


Although I was able to obtain satisfactory results for transmission and reflection (at least their characteristics), the results for permittivity are far from being acceptable. Below I enclosed the neff result for stack with graphene chemical potential set to 1.5 eV. The Re(eps) plot starting at around 7.1 um looks fine, while below this value is “shattered”.

neff visualization


Could you please give me some hints on how to fix this?

Additionally I would like to perform angular analysis of the same structure (also in the same article). I have changed the angle theta property in the source and then tried changing FDTD boundaries from periodic to Bloch or changing plane wave type to BFAST, but to no avail. Either I did not get plane wave any more (for S parameters analysis) or the simulation was diverging. Could you also help me with this issue?

EDIT (27.12.2018):
I guess there is some issue with the s-parameters calculation. I have tried to repeat the same simulation without any graphene layers - only 3 um slab of SiO2 set as a dielectric with constant refractive index = 1.45. The results of neff (both real and imaginary parts) for both Bloch/periodic and BFAST plane waves are below.

Bloch/periodic plane wave


BFAST plane wave


I still did not manage to run this simulation for different angles.

Hi @michal.dudek

Sorry for a late reply. I was on new year break!

Thanks for sharing the screenshots and the new edits. Can you please share your simulation file for the 3um oxide slab and elaborate on how you are calculating effective index i.e. the attached plots?


Dear @bkhanaliloo, I hope you enjoyed your break!

The simulation file with SiO2 slab is here.
The effective index (as well as permittivity) were calculated using the S parameters analysis group from Object Library.

I guess this approach will not be suitable for angular analysis, so I would also be grateful for any hints on how to perform those simulations - especially with BFAST plane wave.

EDIT (07.01.2019):
I think that discontinuities in neff graph are connected with this discussion. Although, I still do not know how to fix this issue.

Hi @michal.dudek

Thanks for the wait and also updating me with your latest findings.

In the analysis group (or the link that you shared with me), n is defined as below:

n = (imag(n1) >= 0)*n1 + (imag(n1) < 0)*n2;

The goal of this calculation is to make sure that imaginary part of effective index is positive so that material does not act as gain material. Thus, the value of n (or effective index neff) depends on the sign of imag(n1). The discontinuity that you are seeing happens when imag(n1) changes the sign. I checked this by choosing n1 and n2 in the results section of s_params analysis group and plotting them:


So there are few things to point out here:

  1. The equation that we are using to calculate neff is not general i.e. you may need to use different equations for different scenarios. This is explained in the references in the provided link. In your case you need to check the paper and see how they have calculated neff specially for graphene that has parallel and perpendicular index.

  2. You need to perform convergence testing and make sure that the results of FDTD simulations converge. If the results are not converging, post processing the data to calculate neff will cause additional issues. For example, in the case where you had a simple lossless rectangle, the imaginary part of n1 and n2 should be zero in the plot above. You can use stackrt to calculate the exact R and T results for the stack of graphene+other material and then use those values in the s_params analysis group.

  3. Sometime the discontinuity can be the result of 2\(\pi\) phase shifts which happens if the wavelength is smaller than the thickness of the material. In this case you need to manually add multiples of 2\(\pi\) phase to n1 and n2 to fix the discontinuity.

Hope this helps.

Hi @bkhanaliloo,

Thank you for your help regarding neff calculations. The discontinuities in my case were coming mostly from calculations of z1 and z2 and then combining them into z. Even z1 and z2 had discontinuities, but I have managed to get rid of them based on the phase unwrapping - the discontinuities were in the places where the phase was wrapped (although with π shifts).

I will try to modify the simulations for plane wave with an angle, and will let you know about the results.

Sure thing. Please feel free to share your results and also the steps that you took to resolve the issue for future references.

Hi @bkhanaliloo,

Regarding the discontinuities in neff calculations - they were originating from the stack structure itself (multiple stacks of SiO2 + graphene). When I have run the simulation for only one layer (SiO2 + graphene), I have obtained good results.

Regarding the angled simulations, I still did not get reasonable results. Could you please check the simulation file and help me with this issue? I know that the s_params analysis group was designed for normal incidence, but the T and R parameters are obtained straight from the T and R monitors, so it should be fine.
In the file I am using 20 stacks of dielectric (n = 1.45) + 2D graphene in the wavelength range from 1 to 10 um (500 points). The angle of incidence is 30 deg for 90 deg polarization. It is already in the Analysis mode.

Hi @michal.dudek

The s-param analysis group in your example is the 2012 version and does not include the angled injection. However, there is an updated 2018 version in the metalens example where includes the angled injections as long as input source is only S or P polarized:

Here is also a good link to learn more about the s_param_adv and its normalization:
Metalens example

Hi @bkhanaliloo,

Thank you for the information about new analysis group. I think it should be also added to the Object Library.
I have tried to run the simulations with new analysis group and there is still an issue with the T and R calculation. I have uploaded new simulation file in Analysis Mode and I would be glad if you could check it. The T and R for theta = 30 deg should look similar to the ones for normal incidence, only blueshifted. In this case there are some “peaks” around 4 um.
The only thing I have changed in the analysis group is to add in a Setup Script these lines for plane wave source:

set("multifrequency beam calculation",1);
set("number of frequency points",500);

Hi @michal.dudek

Thanks for the comment. We plan to add it as we are developing our app gallery. We may also modify the s-param-adv analysis group a bit further.

I think first things would be to make sure that simulation triggers auto shutoff and lower the auto shutoff level.

You can also check your results with stackrt. Please note that for stackrt you will need the graphene effective index (volumetric graphene permittivity model).

It is possible that angled injection excites some of the graphene modes. Do you expect to see resonance modes in graphene around the 4 um range?


Hi @bkhanaliloo

Thank you for your comment. I will try to lower the auto shutoff and repeat the simulation.

The results of R and T obtained for angle 30 degree and 90 deg (TM) polarization are below.

simulation results


You can compare them with results obtained with stackrt for 0 deg (TE) and 90 deg (TM) polarization - see below.

stackrt results


At around 4 um we expect change from elliptical region to hyperbolic dispersion (respectively transmission region and reflectance region) - as seen in the second image.

The T_Gn and R_Gn parameters in the first image are ending at around 3.7 um - at the same wavelength the “theta vs wavelength” goes above 90 deg when the “multifrequency beam calculation” is turned off in the source options. I have seen the post by @fgomez, but did not manage to successfully change the time of the pulse.
Do you have any suggestions regarding this discontinuity both in R and R_Gn curves at around 3.7 um?

Keep in mind that the dielectric material was set to constant refractive index = 1.45 and therefore there is no “dip” coming from silica at around 9 um.

EDIT 13-02-2019:
I have run the simulation for lowered auto shutoff (1e-6) and even extended simulation time to 20000 fs, but it did not change anything in the results.
Then I have removed the mesh override region and changed FDTD mesh settings to “auto non-uniform” with “mesh accuracy” set to 6, “min mesh step” to 0.00005 um and “mesh refinement” to “conformal variant 1”. This change lowered the discontinuity peak at around 3.7 um, but the rest of the R and T curves remained unchanged, as you can see below. So there is still some issue with the angle and wavelength dependency.

new simulation results


After changing “plane wave type” to “BFAST” I did get much better results, as you can see below.

BFAST simulation results


The curves of R and T much better resemble the ones obtained from stackrt, although their peak levels are above 1. Is that an issue with some scaling or normalization?

I am curious about one thing in the “s_params_adv” group - the “S_polarization” result consists of two “polarization” components. Is it therefore possible to obtain both TE and TM polarization in one simulation?

Hi @michal.dudek

Thank you for all the hard work and providing the details about your simulations.

I do not know why I missed it, but for broadband angled injection you should use BFAST. I do not think that the transmission above 1 is due to normalization, however, might be due to early shutoff or poor PML absorption. A good test would be to focus on the wavelengths where transmission is above 1 and run narrowband simulations or even single frequency simulation with Bloch boundary conditions.

You can only have one polarization for the source. However, the analysis group uses gratingpolar script command that returns the results for both S- and P-polarizations.

Hi @bkhanaliloo

Thank you for clarification.

Shutoff level was previously set to 1e-6 and number of PML layers was set to 64 (both for “standard” and “steep angle”). What seems to solve the issue was adding additional mesh override region around the structure with mesh size set to 1 nm.

new R & T results with mesh


The T levels are still above 1, but only slightly. Do you think there may be something else that could improve those results?

Hi @michal.dudek

The results looking great, thanks for all the hard work and reporting back.
I think you can try the convergence testing. Do you see any improvement when you reduce the source bandwidth?


Dear @bkhanaliloo

As you saw previously, the simulations that I performed were for the range 1-10 um. When I tried to change the wavelength range e.g. to 1-9 um or 2-9 um, the simulation started diverging. I did not change anything else besides the wavelength range.
Changing it back to the range 1-10 um made the simulation run normally again.
What may be the cause of this divergence?

Hi @michal.dudek

Most probably this is due to material fitting. Please make sure that “make fit passive” and “improve stability” boxes are checked.

To confirm that divergence is due to material fit, you can enable “specify fit range” and use a wavelength span of 1-10um for it while the source bandwidth is set to 1-9um and see if the simulations converges. If they were still diverging, the issue will be in the boundary condition or stability factor.

I would also try to reduce the source bandwidth to 1-5um. As always, reducing the auto shutoff level or stability factor might also help with the results.

Please keep me updated with your findings.

Hi @bkhanaliloo

Thank you, setting material properties to "make fit passive" = 1 solved the issue with diverging simulation. I did not change the value of "improve stability", as then the fit no longer was in good agreement with material data.

Unfortunately, I have another problem. After upgrading FDTD Solutions to the newest version (v8.21.1882) the simulations for different values of conductivity scaling parameter for graphene material (1 to 6) still produced the same results as simulation with conductivity scaling set to 6. When I checked the properties in Material Explorer, the values were changed to the new one, but it seems that the change was not recognized during simulation. Rolling back to the previous version (v8.21.1854) solves the problem and simulations are correct.

Also the new Visualizer has some issues - e.g. to view several plots at once you have to first click on the ones you want to see, and then you cannot unselect them - previously all plots were shown as default.

Hi @michal.dudek

I ran the graphene_THz_metamaterial.fsp from the example below for “conductivity scaling” of 1 and 4 using 8.21.1882 version:

and the results were different:

Can you please test this example, and upload your simplified simulation file for a review?

There are a few bugs regarding the visualizer that are reported for the fix. Thanks for letting me to know about this issue.


Hi @bkhanaliloo

I heve run the simulations you provided independently on two computers with the FDTD Solutions v. 8.21.1882 and the results were the same in both cases.
You can find the files after simulation here: scaling = 1, scaling = 4.
Is there something I am doing wrong?

Hi @michal.dudek

I think the material fitting causes the issue. After I downloaded the scaling = 1 simulation file, here is what I saw for the material fitting:


After I modified the materiel fit, I guess I got correct results:



Can you please confirm that this solves the problem?