Multipole Decomposition

01/2021: updated to pyGDM v1.1+

New in v1.0.8.

!!CAUTION!! : this is a beta-functionality which has not yet been tested in a sufficient number of different scenarios. The extinction coefficients due to the quadrupole moments seem to be incorrect. While the dipole moments are supposed to work correctly in vacuum, there are still doubts about the results in non-vacuum environment cases.

Here we demonstrate the multipole decomposition of the fields and the extinction spectra of arbitrary shaped nanostructures. The implementation in pyGDM is following reference [1].

In this example we reproduce the case of a silicon nano-cube with side-length 160nm in vacuum, shown in figure 3a of reference [1].

[1] Evlyukhin, A. B., Reinhardt, C. and Chichkov, B. N. Multipole light scattering by nonspherical nanoparticles in the discrete dipole approximation. Phys. Rev. B 84, 235429 (2011)

Simulation setup

from __future__ import print_function, division

## --- load the modules
import numpy as np
import matplotlib.pyplot as plt

from pyGDM2 import structures
from pyGDM2 import materials
from pyGDM2 import fields
from pyGDM2 import core
from pyGDM2 import propagators
from pyGDM2 import visu
from pyGDM2 import tools
from pyGDM2 import linear

Config and run simulation

Si nanocube of 160 nm side-length in vacuum, plane wave illumination.

We discretize with 10 nm to achieve a relatively fast inversion for this demonstration.

## --- simulation initialization ---
mesh = 'cube'
step = 10
geometry = structures.rect_wire(step, L=int(np.round(160/step)),
                                      H=int(np.round(160/step)), mesh=mesh)
material = materials.silicon()
struct = structures.struct(step, geometry, material)

## --- incident field: lin. pol plane wave
field_generator = fields.planewave
wavelengths = np.linspace(530,830,31)
kwargs = dict(theta=0.0, kSign=-1)
efield = fields.efield(field_generator,
               wavelengths=wavelengths, kwargs=kwargs)

## --- environment
dyads = propagators.DyadsQuasistatic123(n1=1, n2=1)

## --- create simulation instance
sim = core.simulation(struct, efield, dyads)

print("N dp={}".format(len(geometry)))

## --- run the main simulation ---
#sim.scatter(method='cupy')   # run on CUDA-GPU using "cupy" (to be installed via pip)

structure initialization - automatic mesh detection: cube
structure initialization - consistency check: 4096/4096 dipoles valid
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/ UserWarning: 3D data. Falling back to XY projection...
  warnings.warn("3D data. Falling back to XY projection...")
N dp=4096
/home/hans/.local/lib/python3.8/site-packages/numba/core/ UserWarning: Numba extension module 'numba_scipy' failed to load due to 'ValueError(No function '__pyx_fuse_0pdtr' found in __pyx_capi__ of 'scipy.special.cython_special')'.
timing for wl=530.00nm - setup: EE 14108.3ms, inv.: 2838.8ms, repropa.: 2324.6ms (1 field configs), tot: 19521.2ms
timing for wl=540.00nm - setup: EE 6399.0ms, inv.: 1469.4ms, repropa.: 1225.1ms (1 field configs), tot: 9329.2ms
timing for wl=550.00nm - setup: EE 3063.9ms, inv.: 1252.5ms, repropa.: 1369.8ms (1 field configs), tot: 5790.7ms
timing for wl=560.00nm - setup: EE 2341.6ms, inv.: 1208.7ms, repropa.: 1412.1ms (1 field configs), tot: 5065.0ms
timing for wl=570.00nm - setup: EE 2279.7ms, inv.: 1333.6ms, repropa.: 1491.0ms (1 field configs), tot: 5202.0ms
timing for wl=580.00nm - setup: EE 8572.1ms, inv.: 1534.9ms, repropa.: 1298.1ms (1 field configs), tot: 11659.8ms
timing for wl=590.00nm - setup: EE 3260.6ms, inv.: 1261.2ms, repropa.: 1472.7ms (1 field configs), tot: 6108.6ms
timing for wl=600.00nm - setup: EE 2227.7ms, inv.: 1271.5ms, repropa.: 1477.8ms (1 field configs), tot: 5088.4ms
timing for wl=610.00nm - setup: EE 2185.8ms, inv.: 1273.1ms, repropa.: 1340.5ms (1 field configs), tot: 5053.6ms
timing for wl=620.00nm - setup: EE 6813.2ms, inv.: 1543.3ms, repropa.: 1304.1ms (1 field configs), tot: 9919.3ms
timing for wl=630.00nm - setup: EE 2684.3ms, inv.: 1261.3ms, repropa.: 1503.6ms (1 field configs), tot: 5563.0ms
timing for wl=640.00nm - setup: EE 2109.7ms, inv.: 1228.1ms, repropa.: 1426.4ms (1 field configs), tot: 4877.2ms
timing for wl=650.00nm - setup: EE 2203.9ms, inv.: 1247.3ms, repropa.: 1342.6ms (1 field configs), tot: 5035.6ms
timing for wl=660.00nm - setup: EE 6548.6ms, inv.: 1506.3ms, repropa.: 1273.7ms (1 field configs), tot: 9581.2ms
timing for wl=670.00nm - setup: EE 4226.0ms, inv.: 1338.4ms, repropa.: 1487.6ms (1 field configs), tot: 7164.0ms
timing for wl=680.00nm - setup: EE 2217.3ms, inv.: 1263.0ms, repropa.: 1494.0ms (1 field configs), tot: 5088.1ms
timing for wl=690.00nm - setup: EE 2305.3ms, inv.: 1264.0ms, repropa.: 1504.6ms (1 field configs), tot: 5180.9ms
timing for wl=700.00nm - setup: EE 3603.9ms, inv.: 1537.4ms, repropa.: 1332.5ms (1 field configs), tot: 6723.3ms
timing for wl=710.00nm - setup: EE 6415.5ms, inv.: 1555.2ms, repropa.: 1277.6ms (1 field configs), tot: 9521.6ms
timing for wl=720.00nm - setup: EE 2250.0ms, inv.: 1243.8ms, repropa.: 1466.3ms (1 field configs), tot: 5071.2ms
timing for wl=730.00nm - setup: EE 2244.0ms, inv.: 1234.8ms, repropa.: 1463.5ms (1 field configs), tot: 5054.4ms
timing for wl=740.00nm - setup: EE 2323.9ms, inv.: 1306.7ms, repropa.: 1317.1ms (1 field configs), tot: 5209.5ms
timing for wl=750.00nm - setup: EE 7525.4ms, inv.: 1597.4ms, repropa.: 1287.4ms (1 field configs), tot: 10686.0ms
timing for wl=760.00nm - setup: EE 3599.0ms, inv.: 1255.7ms, repropa.: 1475.8ms (1 field configs), tot: 6442.8ms
timing for wl=770.00nm - setup: EE 2327.1ms, inv.: 1209.0ms, repropa.: 1381.0ms (1 field configs), tot: 5022.2ms
timing for wl=780.00nm - setup: EE 2135.5ms, inv.: 1205.2ms, repropa.: 1475.0ms (1 field configs), tot: 4930.5ms
timing for wl=790.00nm - setup: EE 6799.7ms, inv.: 1474.0ms, repropa.: 1285.5ms (1 field configs), tot: 9819.1ms
timing for wl=800.00nm - setup: EE 4879.4ms, inv.: 1245.8ms, repropa.: 1482.6ms (1 field configs), tot: 7712.9ms
timing for wl=810.00nm - setup: EE 2229.3ms, inv.: 1257.3ms, repropa.: 1475.2ms (1 field configs), tot: 5066.6ms
timing for wl=820.00nm - setup: EE 2768.3ms, inv.: 1236.3ms, repropa.: 1477.8ms (1 field configs), tot: 5598.3ms
timing for wl=830.00nm - setup: EE 8065.1ms, inv.: 1575.2ms, repropa.: 1210.8ms (1 field configs), tot: 11117.4ms

Extinction spectra

Now we calculate the multipole decomposition of the extinction and compare it to the extinction spectrum obtained from the full discretization. Note that the blue-most resonance around 550nm is mainly electric-octupole (see PRB 84 235429, 2011), hence not contained in our decomposition into dipole and quadrupole moments.

## -- spectra of extinction sections per multipole moment
wl, spec1 = tools.calculate_spectrum(sim, 0, linear.extinct)
ex, sc, ab = spec1.T
wl, spec2 = tools.calculate_spectrum(sim, 0, linear.multipole_decomp_extinct, quadrupoles=True)
ex_p, ex_m, ex_q, ex_mq = spec2.T

plt.plot(wl, ex, label='extinct')
plt.plot(wl, ex_p, label='p')
plt.plot(wl, ex_m, label='m')
plt.plot(wl, ex_q, label='q')
plt.plot(wl, ex_mq, label='mq')
plt.plot(wl, ex_p + ex_m + ex_q + ex_mq, label='multipole sum', dashes=[2,2])

plt.xlabel("wavelength (nm)")
plt.ylabel("extinction cross section (nm^2)")
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/ UserWarning: Multipole decomposition is a new functionality still under testing. Please use with caution.
  warnings.warn("Multipole decomposition is a new functionality still under testing. " +

plot the main contributing dipole moments

We set linear \(X\) polarized, plane wave illumination, so we expect the electric dipole moment to be mainly oriented along \(X\) and the magnetic dipole moment mainly along \(Y\) (incidence along \(Z\)). This can be easily confirmed by the multipole decomposition:

Note: the dipole moments are determined by the arbitrarily chosen incident field amplitude (in pyGDM usually “1”).

## -- spectra of multipole moments
wl, spec_dpdecomp = tools.calculate_spectrum(sim, 0, linear.multipole_decomp, quadrupoles=True)
p = np.array([decomp[0] for decomp  in spec_dpdecomp])
m = np.array([decomp[1] for decomp  in spec_dpdecomp])

plt.plot(wl, np.abs(p.T[0]), label='|p_x|')
plt.plot(wl, np.abs(m.T[1]), label='|m_y|')
plt.xlabel("wavelength (nm)")
plt.ylabel("dipole moment (a.u.)")
/home/hans/.local/lib/python3.8/site-packages/pyGDM2/ UserWarning: Multipole decomposition is a new functionality still under testing. Please use with caution.
  warnings.warn("Multipole decomposition is a new functionality still under testing. " +