Skip to content

Array of spheres on a slab

Spheres with radius 100 nm, permittivity 4 and chirality parameter 0.05 form a square array with a period of 500 nm. They rest on a 10 nm slab with permittivity 3, so the plane of their centers lies 100 nm above the slab. A plane wave hits the structure from below at 17.5 degrees, with k_y = 0.3 k0 and its electric field along x. The S-matrix of the whole stack gives the transmitted and reflected power. Lengths are in nm; the frequency column is k0 / 2π, the inverse vacuum wavelength in 1/nm.

"""A square array of chiral spheres on a thin slab: transmission and reflection."""

import numpy as np

import treams_rs as tr

# Enlarge for a full spectrum: treams uses 100 frequencies.
FREQUENCIES = 5

# The sweep stops short of 1/350 nm^-1, where the diffraction order with
# k_y = 0.3 k0 + 2 pi / 500 starts to propagate and the lattice sum diverges.
k0s = 2 * np.pi * np.linspace(1 / 600, 1 / 350, FREQUENCIES, endpoint=False)
material_slab = 3
thickness = 10
material_sphere = tr.Material(4, 1, 0.05)
lattice = tr.Lattice.square(500)
radius = 100
lmax = 3

# Rows hold (T, R) for each frequency.
power = np.zeros((FREQUENCIES, 2))
for i, k0 in enumerate(k0s):
    kpar = [0, 0.3 * k0]
    sphere = tr.sphere_tmatrix(
        k0=k0, lmax=lmax, radius=radius, material=material_sphere
    )
    spheres = tr.solve_periodic(sphere, lattice=lattice, kpar=kpar)

    # Ports for the diffraction orders kpar + G with |G| <= 0.02 / nm.
    ports = tr.PlaneWavePorts.diffr_orders(kpar, lattice, 0.02)
    slab = tr.slab(k0=k0, basis=ports, thickness=thickness, material=material_slab)
    # The array's S-matrix refers to the plane through the sphere centers,
    # one radius above the slab.
    gap = tr.propagation(k0=k0, basis=ports, distance=radius)
    array = spheres.to_smatrix(ports)
    total = tr.stack([slab, gap, array])

    # A plane wave from below, with k_par = kpar and its electric field along x.
    kz = np.sqrt(k0**2 - kpar[1] ** 2)
    incident = tr.plane_wave(direction=[0, kpar[1], kz], pol=[1, 0, 0], k0=k0)
    power[i] = total.power(incident)

if __name__ == "__main__":
    print("Transmission T and reflection R at oblique incidence")
    print(f"{'k0/(2 pi)':>13}{'T':>13}{'R':>13}")
    for k0, (t, r) in zip(k0s, power, strict=True):
        print(f"{k0 / (2 * np.pi):>13.6g}{t:>13.6g}{r:>13.6g}")

# Source: https://github.com/tfp-photonics/treams/blob/1f5d0d6ebb007288f28bc9e16f6d266e8b55dc39/docs/examples/array_spheres.py
# Adapted: reduced sizes, sweep ends before the diffraction threshold at 1/350,
# plotting removed.
# MIT license, copyright Dominik Beutel; see LICENSE.treams.
import numpy as np
import treams

k0s = 2 * np.pi * np.linspace(1 / 600, 1 / 350, 5, endpoint=False)
material_slab = 3
thickness = 10
material_sphere = treams.Material(4, 1, 0.05)
lattice = treams.Lattice.square(500)
radius = 100
lmax = 3

tr = np.zeros((len(k0s), 2))
for i, k0 in enumerate(k0s):
    kpar = [0, 0.3 * k0]
    spheres = treams.TMatrix.sphere(
        lmax, k0, radius, [material_sphere, 1]
    ).latticeinteraction.solve(lattice, kpar)

    pwb = treams.PlaneWaveBasisByComp.diffr_orders(kpar, lattice, 0.02)
    plw = treams.plane_wave(kpar, [1, 0, 0], k0=k0, basis=pwb, material=1)
    slab = treams.SMatrices.slab(thickness, pwb, k0, [1, material_slab, 1])
    dist = treams.SMatrices.propagation([0, 0, radius], pwb, k0, 1)
    array = treams.SMatrices.from_array(spheres, pwb)
    total = treams.SMatrices.stack([slab, dist, array])
    tr[i, :] = total.tr(plw)

Output

Transmission T and reflection R at oblique incidence
    k0/(2 pi)            T            R
   0.00166667     0.983848    0.0161521
   0.00190476     0.992482   0.00751827
   0.00214286     0.966577    0.0334227
   0.00238095     0.920532    0.0794679
   0.00261905     0.776883     0.223117

Differences from the treams version

  • solve_periodic(sphere, lattice=lattice, kpar=kpar) solves the array once; treams writes TMatrix.sphere(...).latticeinteraction.solve(lattice, kpar).
  • to_smatrix(ports) turns the solved array into an S-matrix; treams writes SMatrices.from_array(spheres, pwb).
  • propagation(distance=radius) takes the distance along z; treams takes the displacement [0, 0, radius], which treams-rs also accepts.
  • plane_wave takes the full propagation direction. treams takes the parallel wave vector together with the plane-wave basis of the ports.
  • Both scripts stop the sweep before 1/350 nm⁻¹. There the diffraction order with k_y = 0.3 k0 + 2π / 500 nm turns from decaying to propagating, and the lattice sum diverges. treams-rs raises a ValueError at that frequency; treams returns a finite number.