simulation of lossy materials

[MahmutRuzi]

Hello, It seems like a nice code , though I have not tried yet. I am wondering if you have tried simulating lossy materials. For example, plasmonic metal nanostructures. I have tried a few RCWA codes, some fails and some require very high Fourier expansion terms. I wonder if you have any insight.

Besides, it would be great if you can add some references on which your code is based.

thanks

[mpasson]

Dear Mahmut,

Thanks for submitting the issue. Simulating lossy materials or metals should be as simple as to insert the complex dielectric function instead of the real one (positive imaginary part means loss, negative imaginary part means gain).

I do not have an extensive experience with lossy material, but I did in the past some calculations (unfortunately not published) with lossy dielectrics and metals to calculate absorption in solar cells, and the results were reasonable. Regarding convergence, I do not expect a big difference between a transparent and lossy dielectrics. Off course the situation is different when dealing with metals, where a much greater number of plane waves is needed to get convergence compared to dielectrics. In metals the problem is mostly related to the fact that the field near a metal surface can break the diffraction limit, and the problem will present itself even with a finite difference method (a much finer grid is usually needed for metal than dielectric).

A good starting point may be the following code snippet, which calculates the fundamental guided mode for a slab of lossy dielectric. If you try it, please be sure to use the most recent version of the code. I may assist better if you describe more in details the problem you need to solve.

# creation of layer
cr = A_FMM.Creator()
cr.slab(12.0+1.0j, 2.0, 2.0, 0.3)
lay = A_FMM.Layer(Nx=50, Ny=0, creator=cr)
lay.transform(ex=0.5)

# coordinate array for plotting
x = np.linspace(-2.0, 2.0, 1001)
z = np.linspace(0.0, 5.0, 501)

# plotting real and imaginary part of dielectric constant
fig, ax = plt.subplots(1,3, figsize=(16, 4))
eps = lay.calculate_epsilon(x=x, z=z)
_=ax[0].contourf(z,x, np.squeeze(eps['eps']).real, levels=41, cmap='Greys')
fig.colorbar(_, ax=ax[0])
ax[0].set_title('Real Epsilon')
_=ax[1].contourf(z,x, np.squeeze(eps['eps']).imag, levels=41, cmap='Greys')
fig.colorbar(_, ax=ax[1])
ax[1].set_title('Imag Epsilon')

# solving for the mode
lay.mode(1.0/1.55)
effective_index = lay.get_index()
print(effective_index[0])

# plotting field
u = lay.create_input({0:1.0})
field = lay.calculate_field(u, x=x, z=z)
ax[2].contourf(z,x, np.squeeze(field['Ey']).real, levels=41)
ax[2].set_title('Fundamental mode field')

# making plots nicer
for a in ax:
    a.set_xlabel('z')
    a.set_ylabel('x')

Regarding some references to papers, I definitely agree with you that they should be there. It is a work in progress.

Best,

Marco

[MahmutRuzi]

Hello Marco, Thank you for the timely and detailed response. You are right, most simulation methods struggle with the simulation of metallic nanostructures. I used meep to do FDTD, but it needs too long, and fine simulation that takes too long. I tried RETICOLO, another RCWA package (Matlab), and it does converge for a specific simulation for Nx=Ny=20. It is fairly quick (~30 mins using i9-11900F @ 2.50GHz). But I wanted to use open-source packages written in Python or Julia, and hope to contribute along the way. So far I tried a few packages but did not get reasonable results or there are not enough details about the codes to be useful. The system I am interested in is plasmonic gold nanobar arrays whose resonance peak appears in the mid IR region. Specifically, I am trying to reproduce Figure 3 and Figure 8 of this published work. H. Altug, et al. [Opt. Express 19(12), 11202–11212 (2011). https://doi.org/10.1364/OE.19.011202

cheers

[mpasson]

Hi Mahmut,

that indeed seems like a tricky problem. The code should in theory be able to handle it, but I have no idea how many Fourier components you would need in the expansion. I added to the documentation a simple example showing the calculation of transmission at normal incidence for a 2D dielectric grating (https://a-fmm.readthedocs.io/en/latest/jupyter/2D_grating.html). This should provide a nice starting point if you want to attempt the simulation of the paper. The main needed modification is the insertion of the frequency dependent dielectric constant of the metal. If let me know if you try, I'll be interested in the outcome, and may provide some help.

A word of advice for speeding up the code. By default, the code can leverage the OMP parallelization of numpy (at least if numpy is installed with pip). This means that most of the operations involving matrices are already parallelized. Usually numpy uses a number of threads equal to the number of processors in the machine (at leas on Linux), but this number can be controlled by setting the environment variable OMP_NUM_THREADS to the desired value. In case of frequency sweeps, it is usually found that the most speedup is obtained by parallelizing most external frequency loop (for example with the multiprocessing module in python). If you do it, make sure that number process launched * OMP_NUM_THREADS < number of processors in the machine. Otherwise you will suffer from very bad performance.

Best, Marco

[MahmutRuzi]

Hello, Dear Marco Thanks for the example code ! I'll try the example and modify it for my needs. I'll let you know how the results look once it is done.

thanks

[MahmutRuzi]

Hello Marco, I had a chance to try out your code. I tried the 2D grating example (https://a-fmm.readthedocs.io/en/latest/jupyter/2D_grating.html). But the results different from what is shown in the example. The transmission is very different and there is an error message: A_FMM/stack.py:354: RuntimeWarning: invalid value encountered in double_scalars np.abs(self.S.S11[j2, j1]) ** 2

Can you please help figure out what the problem is ?

Thanks

[mpasson]

Hi Mahmut,

I can try, but with the amount of information you provided it may not be easy :-).

I can still offer some word of advice:

1. How about the steps before? Are the dielectric profiles of the patterned layer and of the stack displayed correctly? If yes, you could try calculating the modes of each layer to see if they make sense.
2. Have you tried to run directly my notebook? If you fetch the last version of the repository, you'll find the jupyter notebook used for the example in Docs/jupyter/2D_grating.ipynb. If you haven't, I suggest you run directly that.

If you tried both and the problem persist, could you share some info about the system you are using (operating system, python version, libraries versions, ecc), I never had problems across different PCs, but you never know.

Best, Marco

[MahmutRuzi]

Hi Marco, Thanks for the quick response. I tried the jupyter notebook, but encountered similar problem. Searching for the web, it seems general to numpy when deciding a number by very small number close to zero.
I install it using python setup.py install. The error is : python3.10/site-packages/A_FMM-0.1-py3.10.egg/A_FMM/stack.py:354: RuntimeWarning: invalid value encountered in double_scalars I am running it on MacOS Ventura 13.1 (M1 chip). The early steps gives exactly the same result as your example. The only problem is with the transmittance calculations. I'll try to set it up on a ubuntu machine to see if that solves the problem.

Thanks

[mpasson]

OK, that's interesting (and a little problematic).

Unfortunately I do not have access to an M1 Mac for testing. Not sure that is to blame, but if you could just check the results on ubuntu that may help.

Continuing with the Mac system, maybe it is worth to see if the problem is somehow related to the dimensions of the matrices involved. Have you tried other examples? If not, I'll advise the following:

1. Docs/jupyter/1D Scattering Matrix.ipynb: there is no Fourier expansion in this example, only the scattering matrix recursion algorithm. Every matrix involved is 2x2, so it should be well conditioned.
2. Docs/jupyter/PML test.ipynb: Here the Fourier expansion is only 1D, so the matrices involved are much smaller. You can try to play with the truncation order to see if the problem remains.

Also, I would avoid installing the software, as setup.py was never fully tested. Best to just add the repository to your PYTHONPATH environment variable.

If you get any interesting results, please let me know.

[MahmutRuzi]

Thanks for the suggestions. It seems it is due to how I set up the code. I was installing a few packages and somehow mixed it.

Thanks