Skip to content

Cluster of spheres

Four high-index spheres with radii from 75 to 110 nm sit at the corners of a tetrahedron. Multiple scattering couples them: the T-matrix of the cluster includes every interaction and keeps one expansion center per sphere. Expanding it about the origin up to l = 6 gives one global T-matrix, valid outside the sphere of radius 260 nm that encloses the cluster. Rotating the cluster by 90 degrees about y and turning the incident wave with it leaves the cross sections unchanged and rotates the intensity map.

"""Four spheres coupled by multiple scattering: cross sections and near fields."""

import numpy as np

import treams_rs as tr

# Enlarge for a full map: treams uses a 101 x 101 grid.
GRID_POINTS = 7

k0 = 2 * np.pi / 1000
material = tr.Material(16 + 0.5j)
lmax = 3
radii = [110, 90, 80, 75]
positions = (220 / np.sqrt(24)) * np.array(
    [
        [np.sqrt(8), 0, -1],
        [-np.sqrt(2), np.sqrt(6), -1],
        [-np.sqrt(2), -np.sqrt(6), -1],
        [0, 0, 3],
    ]
)

spheres = [
    tr.sphere_tmatrix(k0=k0, lmax=lmax, radius=r, material=material) for r in radii
]
# Cluster.solve couples the spheres and returns one T-matrix with an
# expansion center at every sphere.
tm = tr.Cluster(spheres, positions=positions).solve()

inc = tr.plane_wave(direction=[1, 0, 0], pol="negative_helicity", k0=k0)
xs = tm.cross_sections(inc)

x = np.linspace(-300, 300, GRID_POINTS)
z = np.linspace(-300, 300, GRID_POINTS)
xx, zz = np.meshgrid(x, z, indexing="ij")
grid = np.stack((xx, np.zeros_like(xx), zz), axis=-1)


def intensity_map(tm, inc, radii):
    """Total-field intensity on the grid; NaN where the expansion does not hold."""
    sca = tm.scatter(inc)
    intensity = np.full_like(xx, np.nan)
    valid = tm.valid_points(grid, radii)
    total = inc.efield(grid[valid]) + sca.efield(grid[valid])
    intensity[valid] = 0.5 * np.sum(np.abs(total) ** 2, axis=-1)
    return intensity


intensity = intensity_map(tm, inc, radii)

# One expansion about the origin up to l = 6. Its fields hold only outside
# the sphere of radius 260 that encloses the cluster.
tm_global = tm.in_basis(tr.SphericalBasis.default(6))
xs_global = tm_global.cross_sections(inc)
intensity_global = intensity_map(tm_global, inc, [260])

# Rotating the cluster by 90 degrees about y and illuminating it along -z
# reproduces the first illumination in the rotated frame.
inc_rotate = tr.plane_wave(direction=[0, 0, -1], pol="negative_helicity", k0=k0)
tm_rotate = tm_global.rotate(0, np.pi / 2)
xs_rotate = tm_rotate.cross_sections(inc_rotate)
intensity_rotate = intensity_map(tm_rotate, inc_rotate, [260])

if __name__ == "__main__":
    print("Cross sections (nm^2)")
    headers = ("expansion", "scattering", "extinction", "absorption")
    print("".join(f"{h:>13}" for h in headers))
    for label, result in [("local", xs), ("global", xs_global), ("rotated", xs_rotate)]:
        print(
            f"{label:>13}"
            + "".join(f"{v:>13.6g}" for v in (*result, result.absorption))
        )
    maps = [
        ("local", intensity),
        ("global", intensity_global),
        ("rotated", intensity_rotate),
    ]
    for label, values in maps:
        print(f"\nIntensity at y = 0, {label} expansion (rows: x, columns: z)")
        for row in values:
            print("".join(f"{v:>11.6g}" for v in row))

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

k0 = 2 * np.pi / 1000
materials = [treams.Material(16 + 0.5j), treams.Material()]
lmax = 3
radii = [110, 90, 80, 75]
positions = (220 / np.sqrt(24)) * np.array(
    [
        [np.sqrt(8), 0, -1],
        [-np.sqrt(2), np.sqrt(6), -1],
        [-np.sqrt(2), -np.sqrt(6), -1],
        [0, 0, 3],
    ]
)

spheres = [treams.TMatrix.sphere(lmax, k0, r, materials) for r in radii]
tm = treams.TMatrix.cluster(spheres, positions).interaction.solve()

inc = treams.plane_wave([k0, 0, 0], 0, k0=tm.k0, material=tm.material)
sca = tm @ inc.expand(tm.basis)
xs = tm.xs(inc)

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

intensity = np.zeros_like(xx)
valid = tm.valid_points(grid, radii)
intensity[~valid] = np.nan
intensity[valid] = 0.5 * np.sum(
    np.abs(inc.efield(grid[valid]) + sca.efield(grid[valid])) ** 2, -1
)

tm_global = tm.expand(treams.SphericalWaveBasis.default(6))
sca = tm_global @ inc.expand(tm_global.basis)

xs = tm_global.xs(inc)

intensity_global = np.zeros_like(xx)
valid = tm_global.valid_points(grid, [260])
intensity_global[~valid] = np.nan
intensity_global[valid] = 0.5 * np.sum(
    np.abs(inc.efield(grid[valid]) + sca.efield(grid[valid])) ** 2, -1
)

inc = treams.plane_wave([0, 0, -k0], 0, k0=tm.k0, material=tm.material)
tm_rotate = tm_global.rotate(0, np.pi / 2)
sca = tm_rotate @ inc.expand(tm_rotate.basis)

xs = tm_rotate.xs(inc)

intensity_global = np.zeros_like(xx)
valid = tm_rotate.valid_points(grid, [260])
intensity_global[~valid] = np.nan
intensity_global[valid] = 0.5 * np.sum(
    np.abs(inc.efield(grid[valid]) + sca.efield(grid[valid])) ** 2, -1
)

Output

Cross sections (nm^2)
    expansion   scattering   extinction   absorption
        local       191950       227233      35283.7
       global       191950       227234      35283.8
      rotated       191950       227234      35283.8

Intensity at y = 0, local expansion (rows: x, columns: z)
   0.492519   0.461788   0.451871   0.451982   0.460042   0.485888   0.523742
   0.386665    0.37545   0.396188   0.389244   0.352223    0.34476   0.384914
   0.398542   0.567635    1.36139    1.65666   0.548466   0.356186   0.363744
   0.498961   0.886037    1.57047    2.29615        nan        nan   0.556188
   0.611584     1.5044        nan        nan    1.88694   0.986848   0.575312
   0.652524    1.30309        nan        nan    1.47173   0.779578   0.539411
    0.64797   0.909768    1.09706    1.09561   0.981412   0.735045   0.547142

Intensity at y = 0, global expansion (rows: x, columns: z)
   0.492517   0.461037   0.451883   0.449963   0.455679   0.486385   0.524034
   0.387159   0.372718        nan        nan        nan   0.350069    0.38486
   0.398657        nan        nan        nan        nan        nan   0.365329
   0.492758        nan        nan        nan        nan        nan   0.542634
   0.610437        nan        nan        nan        nan        nan   0.578162
   0.655731    1.30454        nan        nan        nan   0.764384   0.539154
    0.64818   0.905497    1.09644    1.09049   0.982361   0.735165   0.546885

Intensity at y = 0, rotated expansion (rows: x, columns: z)
    0.64818   0.655731   0.610437   0.492758   0.398657   0.387159   0.492517
   0.905497    1.30454        nan        nan        nan   0.372718   0.461037
    1.09644        nan        nan        nan        nan        nan   0.451883
    1.09049        nan        nan        nan        nan        nan   0.449963
   0.982361        nan        nan        nan        nan        nan   0.455679
   0.735165   0.764384        nan        nan        nan   0.350069   0.486385
   0.546885   0.539154   0.578162   0.542634   0.365329    0.38486   0.524034

Differences from the treams version

  • Cluster(spheres, positions=...) holds the uncoupled spheres and solve() couples them; treams writes TMatrix.cluster(...).interaction.solve().
  • cross_sections returns named scattering, extinction and absorption values; treams xs returns the tuple (scattering, extinction).
  • in_basis expands the T-matrix about new centers; treams calls it expand.
  • Each result has its own name (xs_global, intensity_rotate, ...); treams reuses xs, sca and intensity_global.
  • One function, intensity_map, computes all three maps.