Skip to content

Array of spheres from cylindrical waves

Spheres with radius 100 nm, permittivity 4 and chirality parameter 0.05 repeat every 500 nm along z and form a chain. to_cylindrical describes the chain by cylindrical waves, and chains repeated every 500 nm along x make a square array. The array lies in the zx plane, on a 10 nm slab with permittivity 3 that faces the y direction. A plane wave hits the structure from below at 17.5 degrees, with k_x = 0.3 k0 and its electric field along z. This is the structure of the array example, turned so that the chains run along z. The transmission and reflection differ from that example by at most 5e-6, because three axial orders describe each chain. 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 slab, built from cylindrical waves."""

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_x = 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)
pitch = 500
# The array lies in the zx plane: chains of spheres along z repeat along x, and
# the slab and the ports face the y direction.
lattice = tr.Lattice.square(pitch, "zx")
radius = 100
lmax = mmax = 3

# Rows hold (T, R) for each frequency.
power = np.zeros((FREQUENCIES, 2))
for i, k0 in enumerate(k0s):
    # The Bloch vector (k_z, k_x) of the array.
    kpar = [0, 0.3 * k0]

    # A chain of spheres along z, in cylindrical waves.
    sphere = tr.sphere_tmatrix(
        k0=k0, lmax=lmax, radius=radius, material=material_sphere
    )
    spheres = tr.solve_periodic(sphere, lattice=pitch, kpar=kpar[0])
    basis = tr.CylindricalBasis.diffr_orders(kpar[0], mmax, pitch, 0.02)
    spheres_cw = spheres.to_cylindrical(basis)
    # Chains repeated along x make the square array.
    chains = tr.solve_periodic(spheres_cw, lattice=pitch, kpar=kpar[1])

    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 = chains.to_smatrix(ports)
    total = tr.stack([slab, gap, array])

    # A plane wave from below (negative y) with its electric field along z.
    ky = np.sqrt(k0**2 - kpar[1] ** 2)
    incident = tr.plane_wave(direction=[kpar[1], ky, 0], pol=[0, 0, 1], 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_tmatrixc.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)
pitch = 500
lattice = treams.Lattice.square(pitch)
radius = 100
lmax = mmax = 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(pitch, kpar[0])
    cwb = treams.CylindricalWaveBasis.diffr_orders(kpar[0], mmax, pitch, 0.02)
    spheres_cw = treams.TMatrixC.from_array(spheres, cwb)
    chain_cw = spheres_cw.latticeinteraction.solve(pitch, kpar[1])

    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(chain_cw, 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.00751821
   0.00214286     0.966579    0.0334212
   0.00238095     0.920536    0.0794637
   0.00261905     0.776879     0.223121

Differences from the treams version

  • treams keeps the ports in the xy plane and turns the basis internally in SMatrices.from_array. treams-rs needs ports in the zx plane for an array of cylinders, so the script states the turned geometry: Lattice.square(pitch, "zx"), the direction [0.3 k0, k_y, 0] and the field along z. treams gives the parallel wave vector [0, 0.3 k0] in the xy plane and the field along x.
  • solve_periodic solves both the chain and the array of chains; treams writes latticeinteraction.solve(pitch, kpar) twice.
  • to_cylindrical(basis) and to_smatrix(ports) convert the solved arrays; treams writes TMatrixC.from_array and SMatrices.from_array.
  • Both scripts stop the sweep before 1/350 nm⁻¹, where the diffraction order with k_x = 0.3 k0 + 2π / 500 nm starts to propagate and the lattice sum diverges.
  • The results agree with treams to a relative 7e-9. treams sums the row of chains with an automatic Ewald split (the parameter that divides the sum into a real-space and a reciprocal-space part) that leaves errors of about 8e-9. With an explicit split of 0.5 to 1, treams agrees with treams-rs to 1e-14. One explicit split cannot serve the lower frequencies, where the axial orders include decaying waves; see differences.