Internal H-field

Author: Clément Majorel (internal H-field calculation code by C. Majorel)

In this example we calculate the internal magnetic field distribution, trying to reproduce the physics reported by Bakker et al. [1].

First we load pyGDM, construct the geometry and setup the simulation.

[1]: Bakker, R. M. et al. Magnetic and Electric Hotspots with Silicon Nanodimers. Nano Lett. 15, 2137–2142 (2015). (

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

## --------------- Setup
mesh = 'hex'
step = 20.0

R_nm = 90.
H_nm = 160
R = R_nm//step
H = H_nm//step
GAP = 30.

## geometry
geo2 = structures.nanodisc(step, R, H, mesh=mesh)
geo2 = structures.center_struct(geo2)
geo2.T[0] += R_nm+GAP/2.

geo1 = structures.nanodisc(step, R, H, mesh=mesh)
geo1 = structures.center_struct(geo1)
geo1.T[0] -= R_nm+GAP/2.

geo = structures.combine_geometries((geo1, geo2), step)

mat = materials.silicon()
struct = structures.struct(step, geo, mat)

## incident field
wavelengths = np.linspace(450,800, 51)
field_generator = fields.plane_wave
kwargs = dict(theta=[0, 90.], inc_angle=0)     # bottom incidence

efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=kwargs)

## environment: air over glass substrate
dyads = propagators.DyadsQuasistatic123(n1=1.5, n2=1)

## simulation initialization
sim = core.simulation(struct, efield, dyads)

visu.structure(sim, projection='xz')

## --------------- run scatter simulation
sim.scatter(calc_H=True, verbose=True)
structure initialization - automatic mesh detection: hex
structure initialization - consistency check: 1162/1162 dipoles valid

 =============== GDM Simulation Information ===============
precision: <class 'numpy.float32'> / <class 'numpy.complex64'>

 ------ nano-object -------
   Homogeneous object.
   material:             "Silicon, Palik"
   mesh type:            hexagonal compact
   nominal stepsize:     20.0nm
   nr. of meshpoints:    1162

 ----- incident field -----
   field generator: "plane_wave"
   51 wavelengths between 450.0 and 800.0nm
   2 incident field configurations per wavelength

 ------ environment -------
   n3 = constant index material, n=(1+0j)  <-- top
   n2 = constant index material, n=(1+0j)  <-- center layer (height "spacing" = 5000nm)
   n1 = constant index material, n=(1.5+0j)  <-- substrate

 ===== *core.scatter* ======
   NO self-consistent E-fields
   NO self-consistent H-fields
/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...")
/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=450.00nm - setup: EE 5187.5ms, HE 4613.2ms, inv.: 720.5ms, repropa.: 987.8ms (2 field configs), tot: 11527.9ms
timing for wl=457.00nm - setup: EE 2324.4ms, HE 550.6ms, inv.: 78.1ms, repropa.: 371.6ms (2 field configs), tot: 3344.5ms
timing for wl=464.00nm - setup: EE 2298.3ms, HE 552.8ms, inv.: 75.6ms, repropa.: 358.8ms (2 field configs), tot: 3305.8ms
timing for wl=471.00nm - setup: EE 2306.1ms, HE 563.7ms, inv.: 77.3ms, repropa.: 364.7ms (2 field configs), tot: 3332.0ms
timing for wl=478.00nm - setup: EE 2315.3ms, HE 583.7ms, inv.: 77.8ms, repropa.: 374.9ms (2 field configs), tot: 3372.4ms
timing for wl=485.00nm - setup: EE 2246.7ms, HE 583.8ms, inv.: 75.9ms, repropa.: 363.9ms (2 field configs), tot: 3292.0ms
timing for wl=492.00nm - setup: EE 2324.9ms, HE 536.1ms, inv.: 72.9ms, repropa.: 381.7ms (2 field configs), tot: 3335.5ms
timing for wl=499.00nm - setup: EE 2306.8ms, HE 523.3ms, inv.: 73.6ms, repropa.: 364.8ms (2 field configs), tot: 3288.2ms
timing for wl=506.00nm - setup: EE 2276.3ms, HE 525.1ms, inv.: 80.7ms, repropa.: 367.4ms (2 field configs), tot: 3269.4ms
timing for wl=513.00nm - setup: EE 2326.8ms, HE 560.7ms, inv.: 79.0ms, repropa.: 371.3ms (2 field configs), tot: 3358.9ms
timing for wl=520.00nm - setup: EE 2333.5ms, HE 566.7ms, inv.: 77.1ms, repropa.: 369.0ms (2 field configs), tot: 3365.6ms
timing for wl=527.00nm - setup: EE 2318.9ms, HE 569.2ms, inv.: 86.4ms, repropa.: 370.7ms (2 field configs), tot: 3366.4ms
timing for wl=534.00nm - setup: EE 2312.6ms, HE 530.3ms, inv.: 81.4ms, repropa.: 365.8ms (2 field configs), tot: 3310.7ms
timing for wl=541.00nm - setup: EE 2306.0ms, HE 570.8ms, inv.: 78.6ms, repropa.: 364.9ms (2 field configs), tot: 3341.0ms
timing for wl=548.00nm - setup: EE 2322.3ms, HE 524.7ms, inv.: 73.8ms, repropa.: 362.6ms (2 field configs), tot: 3303.9ms
timing for wl=555.00nm - setup: EE 2319.9ms, HE 545.9ms, inv.: 73.3ms, repropa.: 363.9ms (2 field configs), tot: 3323.5ms
timing for wl=562.00nm - setup: EE 2317.1ms, HE 551.1ms, inv.: 73.6ms, repropa.: 368.4ms (2 field configs), tot: 3330.3ms
timing for wl=569.00nm - setup: EE 2308.4ms, HE 542.8ms, inv.: 77.1ms, repropa.: 360.1ms (2 field configs), tot: 3313.7ms
timing for wl=576.00nm - setup: EE 2320.6ms, HE 533.6ms, inv.: 76.4ms, repropa.: 370.1ms (2 field configs), tot: 3322.5ms
timing for wl=583.00nm - setup: EE 2312.0ms, HE 538.5ms, inv.: 75.9ms, repropa.: 362.6ms (2 field configs), tot: 3309.7ms
timing for wl=590.00nm - setup: EE 2306.0ms, HE 529.7ms, inv.: 75.9ms, repropa.: 361.9ms (2 field configs), tot: 3293.4ms
timing for wl=597.00nm - setup: EE 2298.4ms, HE 512.9ms, inv.: 74.5ms, repropa.: 370.4ms (2 field configs), tot: 3275.9ms
timing for wl=604.00nm - setup: EE 2293.9ms, HE 510.4ms, inv.: 77.1ms, repropa.: 370.6ms (2 field configs), tot: 3271.9ms
timing for wl=611.00nm - setup: EE 2316.8ms, HE 519.9ms, inv.: 73.2ms, repropa.: 358.7ms (2 field configs), tot: 3288.7ms
timing for wl=618.00nm - setup: EE 2297.7ms, HE 515.7ms, inv.: 84.3ms, repropa.: 360.1ms (2 field configs), tot: 3277.0ms
timing for wl=625.00nm - setup: EE 2327.3ms, HE 545.9ms, inv.: 82.1ms, repropa.: 372.6ms (2 field configs), tot: 3348.1ms
timing for wl=632.00nm - setup: EE 2291.5ms, HE 522.6ms, inv.: 81.0ms, repropa.: 365.8ms (2 field configs), tot: 3285.6ms
timing for wl=639.00nm - setup: EE 2324.7ms, HE 505.8ms, inv.: 78.2ms, repropa.: 368.1ms (2 field configs), tot: 3297.7ms
timing for wl=646.00nm - setup: EE 2318.3ms, HE 516.4ms, inv.: 77.2ms, repropa.: 374.2ms (2 field configs), tot: 3308.0ms
timing for wl=653.00nm - setup: EE 2302.0ms, HE 521.2ms, inv.: 77.1ms, repropa.: 382.8ms (2 field configs), tot: 3303.7ms
timing for wl=660.00nm - setup: EE 2328.0ms, HE 533.6ms, inv.: 73.5ms, repropa.: 366.1ms (2 field configs), tot: 3321.3ms
timing for wl=667.00nm - setup: EE 2314.8ms, HE 530.8ms, inv.: 72.9ms, repropa.: 377.1ms (2 field configs), tot: 3315.6ms
timing for wl=674.00nm - setup: EE 2283.0ms, HE 542.3ms, inv.: 79.8ms, repropa.: 364.7ms (2 field configs), tot: 3291.0ms
timing for wl=681.00nm - setup: EE 2302.2ms, HE 572.4ms, inv.: 88.5ms, repropa.: 377.7ms (2 field configs), tot: 3364.3ms
timing for wl=688.00nm - setup: EE 2302.5ms, HE 520.7ms, inv.: 74.9ms, repropa.: 372.5ms (2 field configs), tot: 3291.3ms
timing for wl=695.00nm - setup: EE 2312.3ms, HE 510.7ms, inv.: 74.7ms, repropa.: 376.4ms (2 field configs), tot: 3293.9ms
timing for wl=702.00nm - setup: EE 2301.3ms, HE 543.5ms, inv.: 77.7ms, repropa.: 362.5ms (2 field configs), tot: 3307.7ms
timing for wl=709.00nm - setup: EE 2306.6ms, HE 524.5ms, inv.: 74.7ms, repropa.: 366.0ms (2 field configs), tot: 3291.8ms
timing for wl=716.00nm - setup: EE 2279.1ms, HE 586.9ms, inv.: 92.3ms, repropa.: 378.4ms (2 field configs), tot: 3360.3ms
timing for wl=723.00nm - setup: EE 2328.3ms, HE 524.1ms, inv.: 74.1ms, repropa.: 363.1ms (2 field configs), tot: 3311.0ms
timing for wl=730.00nm - setup: EE 2326.1ms, HE 544.2ms, inv.: 84.7ms, repropa.: 388.1ms (2 field configs), tot: 3363.8ms
timing for wl=737.00nm - setup: EE 2310.4ms, HE 537.2ms, inv.: 89.7ms, repropa.: 376.1ms (2 field configs), tot: 3333.8ms
timing for wl=744.00nm - setup: EE 2299.1ms, HE 549.0ms, inv.: 75.5ms, repropa.: 373.9ms (2 field configs), tot: 3318.8ms
timing for wl=751.00nm - setup: EE 2341.4ms, HE 523.4ms, inv.: 73.4ms, repropa.: 377.8ms (2 field configs), tot: 3337.1ms
timing for wl=758.00nm - setup: EE 2294.4ms, HE 516.2ms, inv.: 75.8ms, repropa.: 381.7ms (2 field configs), tot: 3288.9ms
timing for wl=765.00nm - setup: EE 2294.0ms, HE 510.5ms, inv.: 74.9ms, repropa.: 383.3ms (2 field configs), tot: 3283.3ms
timing for wl=772.00nm - setup: EE 2312.3ms, HE 516.4ms, inv.: 77.4ms, repropa.: 373.5ms (2 field configs), tot: 3299.8ms
timing for wl=779.00nm - setup: EE 2322.5ms, HE 523.1ms, inv.: 76.2ms, repropa.: 386.4ms (2 field configs), tot: 3331.5ms
timing for wl=786.00nm - setup: EE 2304.6ms, HE 536.8ms, inv.: 77.8ms, repropa.: 365.3ms (2 field configs), tot: 3304.9ms
timing for wl=793.00nm - setup: EE 2309.6ms, HE 536.2ms, inv.: 90.7ms, repropa.: 381.8ms (2 field configs), tot: 3342.0ms
timing for wl=800.00nm - setup: EE 2316.8ms, HE 520.4ms, inv.: 83.5ms, repropa.: 372.2ms (2 field configs), tot: 3317.8ms

Plot spectra and field-amplitudes at resonance

Note that the internal magnetic field simulation is done by passing the calc_H=True kwarg to scatter.

Now we plot the extinction spectra, as well as the absolute field amplitude of E and H-field inside and around the nano-disc dimer.

We reproduce well the results from Bakker et al.

field_configs = [
    dict(theta=0, wavelength=611),
    dict(theta=0, wavelength=632),
    dict(theta=90, wavelength=527),
    dict(theta=90, wavelength=667)
title_list = [
    'pol. X, ED resonance',
    'pol. X, MD resonance',
    'pol. Y, ED resonance',
    'pol. Y, MD resonance',

r_probe = tools.generate_NF_map_XZ(-225, 225, 51,   -175, 250, 51,     0)

plt.figure(figsize=(12,10), dpi=70)

## --- plot spectra
wl, extsp = tools.calculate_spectrum(sim, 0, linear.extinct)
ex, sc, ab = extsp.T
wl, extmp = tools.calculate_spectrum(sim, 0, linear.multipole_decomp_extinct)
p, m = extmp.T

plt.title("X polarization")
plt.plot(wl, ex, label='tot. ext.')
plt.plot(wl, np.abs(p), label='ED')
plt.plot(wl, np.abs(m), label='MD')
plt.xlabel("wavelength (nm)")
plt.ylabel("extinction section (nm^2)")

wl, extsp = tools.calculate_spectrum(sim, 1, linear.extinct)
ex, sc, ab = extsp.T
wl, extmp = tools.calculate_spectrum(sim, 1, linear.multipole_decomp_extinct)
p, m = extmp.T

plt.title("Y polarization")
plt.plot(wl, ex)
plt.plot(wl, np.abs(p))
plt.plot(wl, np.abs(m))
plt.xlabel("wavelength (nm)")
plt.ylabel("extinction section (nm^2)")

## --- plot field-amplitude maps
for i_plot, field_conf in enumerate(field_configs):
    fidx = tools.get_closest_field_index(sim, field_conf)
    Es,Et, Bs,Bt = linear.nearfield(sim, fidx, r_probe)

    # plt.suptitle(title_list[i_plot])
    im = visu.vectorfield_color(Et, show=0, projection='XZ', fieldComp='abs', interpolation='bicubic', cmap='jet')
    plt.colorbar(im, label='$|E|/|E_0|$')
    visu.structure_contour(sim, projection='XZ', color='w', show=0, dashes=[2,2], lw=1.5, zorder=10)
    plt.xlabel('X (nm)'); plt.ylabel('Z (nm)')

    im = visu.vectorfield_color(Bt, show=0, projection='XZ', fieldComp='abs', interpolation='bicubic', cmap='jet')
    plt.colorbar(im, label='$|H|/|H_0|$')
    visu.structure_contour(sim, projection='XZ', color='w', show=0, dashes=[2,2], lw=1.5, zorder=10)

/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. " +
/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. " +
4 4
4 4
4 4
4 4
4 4
4 4
4 4
4 4