Skip to content

Chain of spheres in cylindrical waves

Two metal spheres with permittivity -16.5 + 1i repeat every 300 nm along z: one of radius 75 nm at (-30, 0, -75) nm and one of radius 65 nm at (30, 0, 75) nm. to_cylindrical turns the solved chain into a cylindrical T-matrix with one axis through each sphere, for the seven axial wavenumbers kz + 2πn / 300 nm with |n| ≤ 3. Each axis collects the waves of its own sphere. A plane wave with kz = 0.005 nm⁻¹ travels almost along x, with equal amplitudes of both helicities. Outside the two cylinders around the axes, the cylindrical waves give the field. Between the spheres, the script expands the waves of all cells in regular dipole waves about each point, as in the chain example. The map shows the real part of E_z.

Outside the cylinders, both descriptions give the same field to 5e-9 at x = ±300 nm and to 0.03 at x = 100 nm, 5 nm outside the cylinder of the 65 nm sphere. There the seven axial orders with |n| ≤ 3 converge slowly; with |n| ≤ 6 the difference drops to 0.004.

"""A chain of sphere pairs along z, described by cylindrical waves: its near field."""

import numpy as np

import treams_rs as tr

# Enlarge for a full map: treams uses a 61 x 31 grid.
GRID_POINTS = (7, 5)

k0 = 2 * np.pi / 700
material = tr.Material(-16.5 + 1j)
lmax = mmax = 3
radii = [75, 65]
positions = [[-30, 0, -75], [30, 0, 75]]
lattice = tr.Lattice(300)
kz = 0.005

spheres = [
    tr.sphere_tmatrix(k0=k0, lmax=lmax, radius=r, material=material) for r in radii
]
chain = tr.solve_periodic(
    tr.Cluster(spheres, positions=positions), lattice=lattice, kpar=kz
)

# Cylindrical waves about one axis through each sphere, for the axial
# wavenumbers kz + 2 pi n / 300 with |n| <= 3.1, the diffraction orders of the
# chain. Each axis collects the waves of its own sphere.
bmax = 3.1 * lattice.reciprocal
basis = tr.CylindricalBasis.diffr_orders(kz, mmax, lattice, bmax, 2, positions)
chain_tmc = chain.to_cylindrical(basis)

inc = tr.plane_wave(
    direction=[np.sqrt(k0**2 - kz**2), 0, kz],
    pol=[np.sqrt(0.5), np.sqrt(0.5)],
    k0=k0,
)

x = np.linspace(-300, 300, GRID_POINTS[0])
z = np.linspace(-150, 150, GRID_POINTS[1])
xx, zz = np.meshgrid(x, z, indexing="ij")
grid = np.stack((xx, np.zeros_like(xx), zz), axis=-1)
sca_cells = chain.scatter(inc)


def lattice_field(points):
    """E_z from the waves of all cells, in regular dipole waves about each point."""
    return np.array(
        [
            inc.efield(r)[2]
            + sca_cells.in_basis(
                tr.SphericalBasis.default(1, positions=[r]), kind="regular"
            ).efield(r)[2]
            for r in points
        ]
    )


# Outside the cylinders that enclose the spheres, the cylindrical waves give
# the field.
ez = np.full_like(xx, np.nan, dtype=complex)
outside = chain_tmc.valid_points(grid, radii)
sca = chain_tmc.scatter(inc)
ez[outside] = (inc.efield(grid[outside]) + sca.efield(grid[outside]))[..., 2]

# Between the spheres, the spherical waves of all cells give the field.
distances = np.linalg.norm(grid[..., None, :] - np.array(positions), axis=-1)
between = ~outside & np.all(distances > np.array(radii), axis=-1)
ez[between] = lattice_field(grid[between])

# Outside the cylinders both descriptions hold: compare them there, row by row.
rows = np.all(outside, axis=1)
difference = np.abs(lattice_field(grid[rows].reshape(-1, 3)) - ez[rows].ravel())
difference = difference.reshape(-1, GRID_POINTS[1]).max(axis=1)

if __name__ == "__main__":
    print("Real part of E_z at y = 0 (rows: x, columns: z)")
    for row in ez:
        print("".join(f"{v.real:>13.6g}" for v in row))
    print()
    print("Largest |E_z| difference of the two descriptions outside the cylinders")
    for xi, value in zip(x[rows], difference, strict=True):
        print(f"x = {xi:>6.0f} nm: {value:.2e}")

# Source: https://github.com/tfp-photonics/treams/blob/1f5d0d6ebb007288f28bc9e16f6d266e8b55dc39/docs/examples/chain_tmatrixc.py
# Adapted: reduced sizes, plotting and the copies of neighboring cells removed.
# MIT license, copyright Dominik Beutel; see LICENSE.treams.
import numpy as np
import treams

k0 = 2 * np.pi / 700
materials = [treams.Material(-16.5 + 1j), treams.Material()]
lmax = mmax = 3
radii = [75, 65]
positions = [[-30, 0, -75], [30, 0, 75]]
lattice = treams.Lattice(300)
kz = 0.005

spheres = [treams.TMatrix.sphere(lmax, k0, r, materials) for r in radii]
chain = treams.TMatrix.cluster(spheres, positions).latticeinteraction.solve(lattice, kz)

bmax = 3.1 * lattice.reciprocal
cwb = treams.CylindricalWaveBasis.diffr_orders(kz, mmax, lattice, bmax, 2, positions)
chain_tmc = treams.TMatrixC.from_array(chain, cwb)

inc = treams.plane_wave(
    [np.sqrt(k0**2 - kz**2), 0, kz],
    [np.sqrt(0.5), np.sqrt(0.5)],
    k0=chain.k0,
    material=chain.material,
)
sca = chain_tmc @ inc.expand(chain_tmc.basis)

x = np.linspace(-300, 300, 7)
y = 0
z = np.linspace(-150, 150, 5)
xx, zz = np.meshgrid(x, z, indexing="ij")
yy = np.full_like(xx, y)
grid = np.stack((xx, yy, zz), axis=-1)

ez = np.zeros_like(xx, complex)
valid = chain_tmc.valid_points(grid, radii)
ez[valid] = (inc.efield(grid[valid]) + sca.efield(grid[valid]))[..., 2]
ez[~valid] = np.nan
sca = chain @ inc.expand(chain.basis)
valid = ~valid & chain.valid_points(grid, radii)
vals = []
for r in grid[valid]:
    swb = treams.SphericalWaveBasis.default(1, positions=[r])
    field = sca.expandlattice(basis=swb).efield(r)
    vals.append(inc.efield(r)[2] + field[2])
ez[valid] = vals

Output

Real part of E_z at y = 0 (rows: x, columns: z)
     -1.13538    -0.950371    -0.633754     -0.24811     0.165893
    -0.943137     -0.63381    -0.207348     0.103191     0.382641
    -0.449631          nan      1.03723      0.32233    -0.178329
      8.13685          nan      1.01914          nan     -2.00549
     0.961256     0.243487    -0.287758    -0.199747     -0.67569
     0.103087    -0.217402    -0.516837    -0.709726    -0.832911
    -0.400668    -0.640795    -0.789931    -0.827765    -0.755714

Largest |E_z| difference of the two descriptions outside the cylinders
x =   -300 nm: 5.11e-09
x =   -200 nm: 1.78e-05
x =    100 nm: 3.11e-02
x =    200 nm: 5.77e-06
x =    300 nm: 1.85e-09

Differences from the treams version

  • solve_periodic(tr.Cluster(spheres, positions=positions), lattice=lattice, kpar=kz) solves the chain once; treams writes TMatrix.cluster(spheres, positions).latticeinteraction.solve(lattice, kz).
  • to_cylindrical(basis) converts the solved chain; treams writes TMatrixC.from_array(chain, cwb).
  • in_basis(basis, kind="regular") expands the waves of all cells about a point; treams writes sca.expandlattice(basis=swb).
  • The periodic response has no valid_points, so the script finds the points outside the spheres with NumPy.
  • The treams-rs script also evaluates the spherical waves outside the cylinders and prints how far the two descriptions differ there.
  • The treams version repeats the map for the neighboring cells with the Bloch phase; both scripts here stop at one cell.