Modeling of particles on a substrate using the dispersion materials of the Drude model

[Anton Utyushev]

Dear all,

I attempt to get scattering spectra for the simple example of a spherical nanoparticle on a metallic substrate, which has diisperesions according to the Drude model in the code itself.

I try to use SMUTHI for this design and I have a problem with convergence.
The problem is that when such a system is simulated, from about 500 to 650 nm unphysical oscillations are observed (see Figure nk_disp.png), but if we set eps constant from any calculated value of eps by Drude formula, a smooth graph is obtained (see Figure nk_const.png).
The question is, how can I fix this error, or maybe I'm doing something wrong that leads to possibly two different values that don't add up. Increasing the l_max value does not lead to any improvement in the solution. I am also attaching the calculation script and graphs. If we need any extras, I can send you whatever you need.

Could someone please help me with my problem, I would be very grateful for any help.

Best regards,
Anton Utyushev

[elek...@gmail.com]

Dear all,

Sorry, the files were not attached, I'll attach them again

Best regards,
Anton Utyushev

среда, 26 октября 2022 г. в 18:06:13 UTC+3, Anton Utyushev:

[Amos Egel]

Dear Anton,

in this mail I'll first suggest a solution for your problem and then add some explanation why the problem occured.

1. Solution

To resolve your issue, I recommend to supply a list of waypoints for a custom Sommerfeld integral contour, for example like this:

waypoints = [0, 0 - 0.1j, 2 - 0.1j, 2, 3]

and then, in the simulation constructor:

simulation = smuthi.simulation.Simulation(layer_system=two_layers,
                                          particle_list=one_sphere,
                                          initial_field=plane_wave,

                                          neff_waypoints=waypoints)

2. Explanation

It seems the issue is related to numerical errors in the Sommerfeld integral that is numerically evaluated when calculating the substrate-mediated self interaction of the particle.

The Sommerfeld integral runs from \kappa = 0 to \kappa = \infinity, where \kappa is the in-plane component of the wave vector in the plane wave expansion of the scattered field. Instead of \kappa, it is more convenient to think of the so-called effective refractive index n_eff, which is defined as \kappa divided by k_0, where k_0 is the vacuum wavenumber. In other words, n_eff = 0 corresponds to a plane wave propagating in z-direction, n_eff = 1 corresponds to a plane wave propagating parallel to the substrate surface and n_eff > 1 corresponds to an evanescent plane wave (in vacuum).

The integrand of the Sommerfeld integral typically shows narrow peaks, which are related to resonances in the layer system, such as waveguide modes or surface plasmons. Due to the peaks, the numerical integration is prone to numerical errors. For that reason, Smuthi automatically evaluates the integral along a deflected contour in the complex plane. The deflected integral contour is chosen such that the vicinity of the waveguide mode and surface plasmon resonances is avoided. That way, a more accurate result can be obtained with the same integral sampling.

However, Smuthi does not exactly evaluate the location of the resonances in the complex n_eff plane. Instead, a rule of thumb is applied: The singularities are assumed to be located in the interval [n_min - 0.1, n_max + 0.2], where n_min is the smallest refractive index in the layer system and n_max is the largest refractive index. So, the deflected integral contour deviates from the real axis only in that interval. In most use cases, this rule of thumbs works well, and the Sommerfeld integral is evaluated with good accuracy.

However, in your use case, there is a resonance that is outside the estimated interval. For example, for a vacuum wavelength of 500nm, the resonance is at n_eff = 1.5, which is larger than n_max + 0.2 = 1.2. Maybe you have a clue about the exact location of the resonance, as your refractive index formula is based on some plasmonic model. Maybe you can calculate, at what where to expect the surface plasmon polariton resonance?

For growing vacuum wavelength, the resonance position wanders towards lower n_eff values, until it finally reaches the interval that is covered by the rule of thumb. I think that is why the unphysical oscillations stopped at some maximal wavelength in your plot.

In our email exchange, you mentioned that you also tried the automatic parameter selection tool, with limited success. Well, the automatic parameter selection tool may decrease the integral sampling, which should eventually yield better results, because even with a narrow peak in the integrand could be met with an integral sampling that is fine enough to sample the peak. However, this  may happen only at very small sampling values.

Now, this is why it makes sense to provide a custom integral contour by defining the waypoints. That way, we extend the interval in which singularities are assumed to n_eff values large enough to cover the resonance in the given use case.

At 500nm, the suggested change of the script leads to an integrand like that:

It is much better behaved and the numerical integral can be expected to have good precision.

Hope this helps, otherwise feel free to ask.

Cheers, Amos

--
Smuthi project repository: https://gitlab.com/AmosEgel/smuthi
Online documentation https://smuthi.readthedocs.io/en/latest/
---
You received this message because you are subscribed to the Google Groups "Smuthi" group.
To unsubscribe from this group and stop receiving emails from it, send an email to smuthi+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/smuthi/602915ae-7e40-4ea1-973b-f8a3f4604e02n%40googlegroups.com.