# Anisotropic medium with spatially varying permittivity

Hello
I want to calculate the modes of a cylindrical hollow waveguide (WG) with inner radius 1 um, and outer radius 1.5 um. The surrounding medium is air (both inner and outer regions). The WG is a composition of at least 10 pairs of alternating metal-dielectric thin (20 nm) hollow cylinders. I could find the modes of this WG with FEEM or MODE. Now I want to employ an effective permittivity tensor, treating the WG material as an homogeneous anisotropic material, so that I can use coarser mesh and save time. Next, I will compare the results from both treatments.

The effective permittivity tensor is given in equation (3) of the article Combined Surface Plasmon and Classical Waveguiding through Metamaterial Fiber Design by Smith et al. This is obtained by geometrically rolling up a planar metal-dielectric bilayer stack (which can be described by parallel and perpendicular effective permittivities, given in equations (1) and (2) in the same article). Therefore, I have an anisotropic material where the permittivity at each location varies as a function of azimuth, phi, which is defined by a piecewise function of x and y coordinates.

I know that Lumerical can easily employ anisotropic tensors if the tensor is diagonal and not dependent on the spatial coordinates. If the permittivity tensor is not diagonal then I need first to make it diagonalised, and introduce the eigenvalues and eigenvectors. Since the permittivities are spatially varying in my case, I also run for loops to calculate the eigenvalues and eigenvectors at every spatial position.

To do these, I modified the example script given in the Lumerical Knowledge Base, on the Matrix Transform page, updating the tensor, and the size of the import object. This script is supposed to calculate the diagonalized effective permittivity tensor and the transformation matrix. Doing so, I suppose I created an import object, with size 3x3x3 um, with my spatially varying permittivity tensor. Since my structure is a hollow cylinder, I added a large ring object (index=1) with inner radius, set to the outer radius of the WG, and a smaller cylinder object (index=1) with radius, set to the inner radius of the WG, so that only the WG region is described by the imported object and the surrounding is overwritten by index 1.

How do these sound? I thought everythingâ€™s cool but the calculated modes with such a setup are not correct. I attached my script and MODE files. Could you please help me understand what is wrong with my setup/script?

matrix_transform.lsf (2.4 KB) WG_eff_perm_660.lms (5.3 MB)

Cheers
Bilge

Thank you for a detailed inquiry and sharing all the supporting files.

Based on the article, your approach looks reasonable to me. Can you please share the MODE and FEEM files where you used stack of metal-dielectric layers to calculate the eigenmodes of a hallow cylinder?

Thanks

Thank you for your reply. I attach the silmulation files and some of the results that Iâ€™ve obtained with them. What I expect to see was to observe similar mode properties to that of the paper at certain wavelengths. For example refer to the Figure 5, they observed a mode which propagates along the cladding at 659 nm, whereas another mode which propagates along the core at 323 nm. There two observations are based on stacked material simulations in COMSOL. On the other hand, refer to the Figure 3c, by employing effective medium, they observed similar modes to those which have been found by stacked material simulations. I would expect the same.

I named the file names systematically. For example,
Eff_perm_350_MODE17_neff_-29961 stands for the mode number 17, calculated by effective material approach at the wavelength 350 nm and the calculated effective index is -2.9961.

It is also confusing that the results from FEEM and MODE simulations for stacked material do not match well ehough, for example at 660 nm. In fact the modes that Iâ€™ve found by MODE at 660 nm does not seem physical at all. And the matching at other wavelengths doesnâ€™t seem satisfactory to me.

I have one more questions. How exactly should I read the TE polarization fraction and waveguide TE/TM fraction results? I suppose that for example if the TE polarization fraction is 0, then the mode is totally TE polarized. And for example lets say waveguide TE/TM fraction is 8/99, what does that mean exactly? How do I interpret these together? Is it possible to calculate those with FEEM?

Thanks and I am looking forward to hear from you.
Cheers
Bilge.

Figures:
Eff_perm_310_MODE01_neff_14458

Eff_perm_350_MODE01_neff_-33168

Eff_perm_350_MODE17_neff_-29961

Eff_perm_660_MODE01_neff_85581

Stack_310_FEEM_MODE01_neff_14453

Stack_310_MODE_MODE01_neff_14407

Stack_350_FEEM_MODE17_neff_09912

Stack_350_MODE_MODE09_neff_09961

Stack_660_FEEM_MODE01_neff_18758

Stack_660_MODE_MODE01_neff_19104

Effective_MODE_660_layout.lms (5.3 MB) matrix_transform.lsf (2.4 KB) Stacked_FEEM_660_layout.ldev (5.8 MB) Stacked_MODE_660_layout.lms (307.5 KB)

1 Like

Thanks again for sharing all the files and a detailed inquiry.

I was expecting that FEEM results converge better than FDE results for curved geometries. FEEM uses finite element method while FDE uses rectangular meshes.

As you have noticed, when you use curved geometries in FDE, you create hot spots and this is why the FDE plots have tiny hot spots with no obvious field profile. One approach would be to modify the color bar limits or use log scale for better visualization.

You can improve the FDE results by refining the mesh. In your case this requires a lot of memory. Using symmetry could also reduce the simulation time and memory requirement if you know the symmetry of more of interest.

Regarding the second inquiry: Please visit the link below for the equations used to calculate these quantities:

If the TE polarization fraction is 0 (100), then the mode is TM (TE) polarized.

I will keep working on your files to see if I can find a way to calculate the modes in FDE, and how we can make the effective index approach to work.

Thanks

Hi @bkhanaliloo
Thank you very much for your detailed explanation.
I think I can rely on my FEEM solutions, rather than the FDE solutions, as the FEEM meshing is much more useful for curved geometries. I donâ€™t think for now it is necessary to put more efforts to obtain reliable results by FDE solutions by for example improving the spatial resolution, as I am able to to do so with FEEM. Iâ€™ve already tried finer mesh, and it indeed required too much computational power.
The crucial thing to me is to get solutions with effective medium approach. So, I will be looking forward to hear from you more about this.
Thanks!
Bilge.

We can use the results of FEEM solver to get an idea about the field profile and also use the value of neff in FDE calculations to help the software to locate the modes. However, unfortunately we will not be able to proceed with FEEM for effective material approach as we are lacking the grid attribute (and I could not find any other workaround).

For the effective material approach in FDE, can you please use higher values for matrix transformation?

Currently you have only 20 points for the x and y directions, while the change in permittivity might require higher number of values to capture the spatial variation. I expect that this will be similar to refining the mesh and we might be lucky to avoid the hot spots (local cavities).

Please keep me updated. We have a long weekend here and will be back on Tuesday.

Dear @bkhanaliloo
I hope you had a nice holiday.
I am not sure I understood what you meant correctly. I was thinking that the graphical rendering options were just about the view of the structures. Do they even have anything to do with the simulations? I tried 100 for each of them but nothing has changed.
Thanks
Bilge.

Thanks!

You are right., the â€śMatrix transformâ€ť is for graphical rendering. To be more precise, if you have a lower resolution for the positions where effective index values and grid attributes are defined, increasing graphical rendering will reveal it as there will be no added matrix transform arrows to the GUI.

With the settings that you have, the x and y resolutions are 25nm while the thickness of each stack is 20nm. If we want to have a grid resolution of 5nm, you need more points for x and y axes (121*5):

# define positions where the object refractive index values and grid attribute will be defined
# a cube where x, y, and z span from -1.5-1.5 um
x = linspace(-1.5e-6,1.5e-6,121*5);
y = linspace(-1.5e-6,1.5e-6,121*5);
z = linspace(-1.5e-6,1.5e-6,3);


I would start with twice your current resolution (121*2) and see how the results get improved. The script for 5 times resolution might be quite slow and I recommend to increase it slowly while running the simulations to calculate the modes. I hope that this will improve the hot spot effects.