PyStratify

Light scattering and emission by multilayered films, concentric cylinders and concentric spheres in Python, under plane waves, electric and magnetic point dipoles and focused (Cassegrain) beams. One modal-response core serves the three geometries: the transfer-matrix method for films (after Byrnes's tmm), its cylindrical analogue for cylinders, and for spheres the recursive transfer-matrix method of Moroz (2005) as used by STRATIFY (Rasskazov, Carney & Moroz, OSA Continuum 3, 2290 (2020)). Spheres go further: absorbing, magnetic, gain and chiral (Pasteur) shells, 2D materials and Feibelman d-parameters on the interfaces, and electric, magnetic, chiral (p + m) and electric-quadrupole emitters inside or outside the particle, with angle-resolved far fields and decay rates. The numerics are rebuilt so that every quantity stays finite and accurate to the orders it needs, and the defects found in the MATLAB code are fixed (AUDIT.md).

0.10.5: cylinder emitters at a mode's cutoff

Next to a light line the reflected spectral density of a cylinder source loses digits: its m = ±1 TE/TM-coupled terms carry 1/q factors that cancel, ~1e-4 ε/(1 − b)² (1e-5 at 1e-7 from the line, noise at 1e-10). A guided mode just above cutoff (HE₁ₙ modes approach it exponentially) puts a sharp feature there, and the quadrature chased the noise until its budget ran out; with scipy 1.18 even a plain fibre did so at tolerance 1e-7. Within 1e-3 of the host's and the source medium's light lines the real part now comes from the power balance at that axial wavenumber, escape + absorbed − direct (direct = ¾(1 + b²) transverse, 3/2(1 − b²) axial), which equals the reflected density to ~1e-14 wherever that is accurate; elsewhere the balance stays an independent check. Converged results are unchanged.

0.10.4: every guided mode of a cylinder

A point source beside or inside a lossless cylinder could lose the power it couples into some of the cylinder's guided modes. The mode search minimised the dispersion determinant on a grid and kept a minimum only if the minimiser's last value was below a fixed bar; at a root |det| is V-shaped, Brent's parabolic steps stalled just above it, and real modes were dropped while the solve reported convergence. Over 60 random lossless cylinders (electric and magnetic sources) the totals moved by more than 1e-3 in 28 of 97 comparable cases (up to 9.5%), and 11 that raised an error now solve. Roots are now bracketed and judged at the root, against the determinant's own scale; the light line, where modes cut off, is sampled geometrically and never evaluated within 1e-12 of; the residue's error estimate measures the value used. For a step-index fibre the modes found are exactly those of its closed-form dispersion relation, and the rates are smooth through a mode's cutoff. Found by a 3D FDTD comparison with Meep (benchmarks/FDTD_MEEP.md). Near fields are now exact at a sphere's centre (below).

0.10.3: guided power from a source outside the stack

A point source in an exterior half-space of a film stack (above or below it) lost the power it couples into the stack's guided modes: the transverse-resonance search ran at the source plane, where the source side has no interface, and found no poles. Total rates near lossless guides were low (by 20%, 30 nm above 200 nm of TiO2 on glass) while the solve reported convergence; lossy stacks, and sources inside a finite layer, were unaffected. The modes are now found inside the stack. smuthi and the lossless limit of an absorbing guide confirm it (benchmarks/EMITTER_REFERENCES.md).

Focused illumination (0.10.0)

A fourth source for all three geometries: a beam focused by an aplanatic objective whose pupil is filled uniformly, through an annulus (a Cassegrain, obscuration > 0) or a disc, as in FT-IR microspectroscopy (Davis, Carney & Bhargava, Anal. Chem. 82, 3474 (2010) and 83, 525 (2011); van Dijk et al., Appl. Spectrosc. 67, 546 (2013)). Powers are fractions of the beam's power; the detected one passes a condenser's annulus (transmission) or the objective's (film transflection), and the apparent absorbance is −log₁₀ of its ratio to a reference setup's.

import numpy as np
import pystratify as ps

beam = ps.FocusedBeam(0.5, 0.23)                       # NA, obscuration NA (n sin θ in the host)
wavenumber = np.linspace(1000, 4000, 301)              # cm⁻¹
wavelength = 1e4 / wavenumber                          # µm
fibre = ps.focused_cylinders([5.0], [1.49 + 0.01j, 1.0], wavelength, beam)
fibre.apparent_absorbance, fibre["extinction"], fibre.by_polarization["y"]["detected"]
film = ps.focused_films([1, n_polymer, 1.43, 1], [np.inf, 2.0, 400.0, np.inf], wavelength, beam,
                        coherence=["i", "c", "i", "i"],
                        reference=([1, 1.43, 1], [np.inf, 400.0, np.inf], ["i", "i", "i"]))
sphere = ps.focused_spheres([5.0], [n_pmma, 1.0], wavelength,
                            [ps.FocusedBeam(0.5, 0.23, offset=(x, 0, 0)) for x in (0, 2, 4)])  # a line scan
field = ps.focused_field_spheres([5.0], [n_pmma[0], 1.0], wavelength[0], beam, points)    # |E|² = 1 at the focus

FocusedBeam(na, obscuration, collection_na, collection_obscuration, mode, illumination, field_stop, model, polarization, offset): the condenser defaults to the objective; illumination="kohler" makes every pupil direction an independent plane wave filling a field stop of diameter field_stop (a cylinder through it is lit over that chord); model="scalar" is the corrected scalar theory of the coherent sphere and cylinder spot; polarization is the pupil's ("x", "y" or "unpolarized"; a cylinder's axis is y); offset moves the focus (coherent only). solve_problem(Problem(..., source=beam)) takes one wavelength, with reference=Problem(...) for films and points= for "field"; the focused_* functions batch wavelengths and take a list of beams (an NA sweep or a focus scan), sharing every Mie, oblique-cylinder and planar solve between them.

• Films: the cone average, uniform in sin²θ, s and p equally, of plane-wave R and T and each layer's absorptance, by adaptive quadrature split at the condenser's edges and every lossless layer's critical angle, with the plane-wave solver the coherence column calls for (planar.layered_response). A 400-µm coherent window is resolved fringe by fringe. Coherent and Köhler illumination coincide.
• Spheres: the beam's vector multipole coefficients (m = ±1, a 1-D pupil integral, when the focus lies on the axis through the centre; every m, by FFT in the azimuth, otherwise), the T-matrix, P_ext = −4 Re Σ T|c|², P_sca = 4 Σ |Tc|² and the collected power over the condenser in closed form per order. Köhler uses the cross sections over the cone and the scattering landing in the condenser (Legendre addition theorem).
• Cylinders (axis y, beam z): each pupil row s_y is a 2-D problem; one oblique solve per row serves every transverse angle, the rows add in power, and the s_y integral is adaptive with an error estimate, split where rows cross an annulus edge. The scattered light is collected over the condenser at its true angle, asin(s_x/√(1−s_y²)) — never masked to the illumination.
• Scalar model: U and ∂U/∂r continuous; the sphere's order-l coefficient is the TE one plus a monopole, the cylinder's the axial-electric one at normal incidence with indices √(n² − (n_host s_y)²); the pupil weighting is the vector model's, so vector − scalar isolates polarization.
• Near fields (coherent): spheres from the beam's multipoles and the shells' radial functions, cylinders as the coherent sum over rows, films as every pupil direction's coh_tmm field summed over the azimuth in closed form (Bessel J₀, J₁, J₂); the beam alone in the host is the Richards–Wolf integral.

The ambient (films) or host must be lossless and NA < n there. pytest tests/test_focused.py checks every semi-analytic form against a plane-wave superposition on a dense pupil grid through amplitude_matrix, the oblique solve_cylinder, coh_tmm/position_resolved, near_field and CylinderSolution.field; energy balance (lossless: absorption 0, per row for cylinders); NA → 0 against the plane-wave cross section and cross width; an index-matched sphere and an empty stack invisible; the aplanatic film path factor ⟨1/cos θ⟩ = 1.0357 (n = 1.5, NA 0.23–0.5); Köhler against cone averages of plane-wave results; the scalar coefficients against closed forms; and the film near field's lateral integral against the cone-averaged plane-wave intensity (Parseval).

Unified films, cylinders and spheres (0.9.0)

solve_problem(Problem(...)) chooses the geometry's modal basis and uses the owned response core. NumPy/SciPy supply numerical primitives; external optical engines are benchmarks. The development contract accepts passive, isotropic, nonmagnetic media, plane waves, and electric/magnetic 3D point dipoles in lossless source and exterior regions. Cylinder point sources integrate the continuous axial spectrum.

import numpy as np
import pystratify as ps

film = ps.Problem("films", [np.inf, 0.2, np.inf], [1, 1.5, 1], 0.6,
                  ps.PointDipole(0.1, "electric", layer=1,
                                 intrinsic_quantum_yield=0.8))
rates = ps.solve_problem(film)
assert rates["diagnostics"]["converged"]

cylinder = ps.Problem("cylinders", [0.05, 0.1], [1.5, 2, 1], 0.6,
                      ps.PlaneWave(angle=np.pi / 3))
widths = ps.solve_problem(cylinder)

Lengths share one caller-selected unit; angles are radians. Film thickness arrays include infinite exteriors, and source positions are depth within the specified region (negative in the incident region). Radial arrays are increasing interface radii with the exterior index last; source positions are absolute radii. Cylinder plane-wave angles are from +z. Widths use total incident irradiance and have length units; sphere cross sections have area units.

Rates include canonical orientations and their isotropic average, total, escape, guided, absorbed, efficiencies and diagnostics. The reference is the same dipole in its own homogeneous medium. lifetime_factor = (1-q0) + q0*total is the inverse lifetime ratio; its reciprocal is the lifetime ratio. Material absorption is distinct from intrinsic nonradiative loss. Check convergence before using results. Numerical budget/precision failures raise exceptions; the HTTP adapter masks those samples.

Lossy cylinder spectra guard narrow intervals next to lossless-medium light lines, where the vector basis loses floating-point precision. The endpoint continuation uncertainty is reported as diagnostics.grazing_error and included in the integration error; it must meet the requested tolerance alongside independent channel balance and angular-tail checks. Lossless guided-mode integration is unchanged. Scaled Hankel functions use direct evaluation with logarithmic high-order overflow recovery, and pattern cuts at the same polar angle share their mode coefficients.

Patterns use outputs=("rates", "pattern"), theta=..., phi=... with broadcast angle arrays in radians, and return normalized power per solid angle. Film theta zero points outward into the incident medium; spheres/cylinders use the +z polar axis. Cylinder waves interpret their theta pattern samples as azimuth. The sphere API below is part of the stable 1.x surface, except shell_average, spectral_density, casimir_polder and green_dyadic, which are provisional, and the normalized-route functions, which are legacy (docs/API.md).

response.py owns scalar, logarithmic and coupled block composition. Geometry bases, source spectra and reconstruction are separate modules. The migrated planar routines retain the extended API and upstream MIT notice in src/pystratify/PLANAR-LICENSE.txt. Coherent propagation has no artificial attenuation clipping. Texture emission, line sources, finite cylinders and absorbing-source regularization are deferred. Run pytest tests/test_unified.py for analytic Fresnel/Bessel limits, independent flux/Green closure, guided channels, patterns, ED/MD, metal-cylinder endpoint error budgets, mode reuse, high-order special functions and treams cylinder benchmarks.

No optical constants are shipped: pass n + ik from a database such as refractiveindex.info. The Drude fits in DRUDE exist only for what tables cannot supply, the free-electron damping.

Install

pip install pystratify

The development version: pip install git+https://github.com/multilayerapp/pystratify.git.

Python ≥ 3.10, NumPy ≥ 1.24, SciPy ≥ 1.10.

Use

Which names are the stable 1.x surface and which are provisional, legacy or internal: docs/API.md (ruled 2026-10-10; binding from 1.0).

Lengths are unit-agnostic (give radii and wavelengths in the same unit); the Drude helpers and graphene_conductivity take nanometres. Shells are indexed from 0 (the core) to N (the host).

import numpy as np
import pystratify as ps

wavelength = np.linspace(500, 900, 401)                              # nm
radii = [50.0, 55.0]                                                 # SiO2 core, Au shell
n = np.stack([np.full(401, 1.45), n_gold, np.full(401, 1.33)], 1)   # (W, N+1), host last
sol = ps.solve(radii, n, wavelength)                                 # mu defaults to 1

cs = ps.cross_sections(sol)                 # cs.q_ext, cs.q_sca, cs.q_abs: arrays over wavelength
cs.ext_by_order[:, ps.TM, 0]                # electric dipole contribution
field = ps.near_field(sol, x, y, z, wavelength_index=200)  # fields at the very surface:
                                            # ps.solve(..., regime="near") converges them
stored = ps.shell_energy(sol, wavelength_index=200)

rates = ps.decay_rates(radii, n[200], wavelength[200], r=[60.0, 75.0])
rates.radiative, rates.nonradiative, rates.total, rates.balance_error, rates.quantum_yield()

theta, phi = np.meshgrid(np.linspace(0, np.pi, 181), np.linspace(0, 2 * np.pi, 361), indexing="ij")
pattern = ps.scattering_pattern(sol, theta, phi, polarization=(1, 1j))  # helicity +1 along +z
pattern.directivity, pattern.differential_cross_section, pattern.stokes()   # (W, 181, 361)
emitter = ps.dipole_far_field(radii, n[200], wavelength[200], position=[0, 20, 60],
                              moment=[1, 0, 0], theta=theta, phi=phi)
emitter.directivity, emitter.power, emitter.circular_polarization          # power = P_rad / P_0

kappa = [0, 0, 1e-3, 0]                     # Pasteur chirality per shell, host achiral
chiral = ps.solve_chiral([50.0, 55.0, 60.0], [1.45, n_gold[200], 1.5, 1.33], kappa, wavelength[200])
cd = ps.helicity_cross_sections(chiral)    # .ext[:, 0 or 1] for helicity +1 / -1, .cd_ext, .g_ext
ps.scattering_pattern(chiral, theta, phi)  # far fields and dipole_far_field(..., kappa=...) work alike

# a chiral molecule (m parallel to p, a quarter period out of phase) inside the chiral shell,
# randomly oriented; and a valley exciton of a 2D semiconductor, a circular in-plane dipole
molecule = ps.dipole_far_field([50.0, 55.0, 60.0], [1.45, n_gold[200], 1.5, 1.33], wavelength[200],
                               [0, 0, 57.5], [0, 0, 1], theta, phi, kappa=kappa,
                               magnetic_moment=[0, 0, 0.01j], orientation="isotropic")
molecule.dissymmetry, molecule.helicity_power, molecule.coherency   # g_lum, (P_+, P_-), <F F^H>
exciton = ps.dipole_far_field(radii, n[200], wavelength[200], [0, 0, 57.0], [1, 1j, 0], theta, phi)

# decay of a randomly oriented chiral molecule at two distances from a chiral sphere; the enantiomer
# has m -> -m.  (A lossy kappa needs magnetic loss too: passivity is eps'' mu'' >= kappa''^2.)
rates = ps.emission_rates([50.0], [1.6 + 0.02j, 1.0], 600.0, [[0, 0, 55.0], [0, 0, 70.0]], [0, 0, 1],
                          magnetic_moment=[0, 0, -0.1j], kappa=[0.03 + 0.003j, 0], mu=[1 + 1e-3j, 1],
                          orientation="isotropic")
rates.total, rates.radiative_helicity, rates.absorption, rates.balance_error, rates.dissymmetry

# an electric-quadrupole emitter 5 nm from gold: rates and far field
quad = ps.emission_rates([40.0], [n_gold[200], 1.33], wavelength[200], [0, 0, 45.0], [0, 0, 0],
                         quadrupole=np.diag([-1.0, -1.0, 2.0]))

# 2D materials on interface 0: a TMD monolayer as a sheet, graphene from its conductivity
tmd = ps.Sheet.from_film(eps_ws2, 0.618, wavelength, eps_background=1.0, eps_normal=6.5)
sol = ps.solve([78.0], [3.9, 1.0], wavelength, sheets={0: tmd})
sigma = ps.graphene_conductivity(wavelength_ir, fermi_energy=0.5, damping=0.005)  # sigma Z0; nm, eV
sol = ps.solve([25.0], [1.45, 1.0], wavelength_ir, sheets={0: sigma})

Far-field conventions are Bohren & Huffman's: E_sca = e^{ikr}/(−ikr) X with (X_∥, X_⊥) = [[S2, S3], [S4, S1]] (E_∥, E_⊥) and dσ/dΩ = |X|²/k². Helicity +1 is (x̂ + iŷ)/√2 for a wave along +z (left-circular in the optics convention). Chiral media follow D = εE + iκH, B = μH − iκE (Gaussian units, e^{−iωt}; the convention of treams), so helicity ±1 has index n ± κ.

Emitters follow Klimov, Guzatov & Ducloy (EPL 97, 47004 (2012)) and Guzatov & Klimov (New J. Phys. 14, 123009 (2012)): a chiral emitter is an electric dipole p and a magnetic dipole m radiating coherently, Gaussian units (in SI pass m/c), enantiomers differing in the sign of Im(p·m*); their molecule ⟨e|d|g⟩ = d₀, ⟨e|m|g⟩ = −i m₀ is (p, m) = (d₀, −i m₀). With both moments the far field is normalised to the electric dipole's, so m enters as (n_h/μ_h) m, and in a host of index n a free emitter has g_lum = 2(P₊ − P₋)/(P₊ + P₋) = 4n Im(p·m*)/(|p|² + n²|m|²). m is the dual of p (the magnetic current −iωm of the Maxwell equations), as in decay_rates; magnetic_convention="current" takes a current-loop moment m_A instead, which in a layer with μ_d, κ_d acts as the dual moment μ_d m_A together with the electric dipole iκ_d m_A. An electric quadrupole is Jackson's Q_ij = ∫(3x_ix_j − r²δ_ij)ρ dV: it couples through (1/6)Q:∇E and alone in the host radiates (k²/120)Σ|Q_ij|² in units of the power of |p| = 1. Rates are normalised to the same source in the unbounded host (normalization="layer": in its own layer's medium). A lossy chirality (circular dichroism) is passive only with magnetic loss, ε″μ″ ≥ κ″²; emission_rates warns otherwise, since one handedness then sees gain. A valley exciton of a monolayer semiconductor is a circular in-plane dipole (x̂ ± iŷ)/√2 about the layer normal (Gong et al., Science 359, 443 (2018)); dark excitons are out-of-plane. For an ensemble of incoherent emitters uniformly covering a sphere of radius r₀ (a conformal monolayer, a molecular shell) with orientation fixed relative to the local normal, the emission is isotropic and its dissymmetry is that of one emitter at r₀ẑ, so power and helicity_power of a single call give the ensemble.

A 2D material on interface j enters as generalised sheet transition conditions with the fields averaged over its two sides (Kuester et al., IEEE TAP 51, 2641 (2003)): r̂ × (H⁺ − H⁻) = σ⟨E_t⟩ and E_t⁺ − E_t⁻ = −ζ∇t⟨D_n⟩, with σ = σ_s Z₀ (SI) = 4πσ_s/c the dimensionless sheet conductivity (graphene's universal value is πα) and ζ a length for the out-of-plane response. Sheet.from_film gives the first-order limit of a film of thickness d, σ = −ik₀d(ε∥ − ε_b), ζ = d(1/ε_b − 1/ε_⊥), which absorbs (c/8π)[Re σ|⟨E_t⟩|² + k₀ Im ζ|⟨D_n⟩|²] per area. graphene_conductivity is the local RPA conductivity (Falkovsky; Hanson), intraband at any temperature, interband for k_BT ≪ E_F.

examples/ has runnable scripts (multipole spectra, free-path correction, scattering pattern, scattering and emission directivity, near-field maps, energy density, stored energy, electric and magnetic dipole decay, circular dichroism of a chiral shell, circularly polarised emission of chiral molecules and valley excitons beside spheres, enantioselective decay beside chiral particles, quadrupole versus dipole emitters near gold, a TMD monolayer on a silicon sphere and graphene-coated spheres); with matplotlib installed each saves a PNG beside itself. crosscheck/ recomputes published figures from their stated parameters for comparison with the papers (Guzatov & Klimov 2012, Figs. 3–6).

Numerics

The transfer matrices of the method are products of Riccati–Bessel functions ψ_l, ξ_l at different arguments. Past the size parameter they over- and underflow, and inside thin metal shells their differences cancel long before that: in the 3-nm gold shell of OSAC Fig. 2(b) a direct product violates the interface conditions by ~50 % at l = 14, while an emitter a few nm from a metal needs l ≈ 500–1000. PyStratify keeps the physics and changes the representation.

• Special functions as complex logarithms (riccati.py). Orders up to |x| + 32 come from the exponentially scaled AMOS routines, which are exact near the zeros of ψ; higher orders from the stable recurrences, downward for ψ and upward for ξ, so the cost does not grow with the truncation order. In gain media AMOS's hankel1e returns spurious zeros and yve is mis-scaled; there the unscaled hankel1 is used. Off the real axis ψ has no zeros, so a zero from AMOS is underflow (scaled jve underflows below n = |x| once Im x ≳ 1000) and the recurrence takes over. Accuracy is ≤ 1e-10 against extended-precision mpmath for orders 0–2500, |x| from 1e-6 to 3000, Im x up to 700 and gain, and 1e-12 at 1500 + 1500i; spheres with |mx| = 3·10⁴ and Im(mx) = 8·10⁴ match BHMIE to 2e-12.
• Scaled two-sided recursion (solver.py). The regular solution is swept outwards through ρ = Rξ/ψ and the outgoing solution inwards through σ = Sψ/ξ, both carried as logarithms, so a ratio that falls to e⁻⁷⁰⁰ across a thick shell keeps full precision instead of turning subnormal. At high order ψ'/ψ ≈ (l+1)/x on both sides of an interface, so the mismatch that sets ρ is a small difference of large numbers; it is never formed by subtraction but assembled from the material contrast and the small ratios ψ_{l+1}/ψ_l, ξ_{l−1}/ξ_l. Each amplitude uses the better-conditioned of the value and derivative matching equations. The approach follows Yang (2003), Peña & Pal (2009), Ladutenko et al. (2017), Majic & Le Ru (2020) and Zhang (2025).
• Decay rates from the Green's function. The radiative rate, the Ohmic loss and the total rate (LDOS) are three independent sums, and DecayRates.balance_error checks total = radiative + nonradiative. The loss includes magnetic loss, Im(μ)|H|², as well as Im(ε)|E|². Sums are truncated from the geometry and accepted only when a remainder estimate meets tol; an unconverged position is flagged, never silently truncated.
• Dipole far fields from reciprocity. The far-field amplitude of a source (p, m) at r₀ in direction n̂ with polarisation ê is p·E(r₀) − m·H(r₀) for the plane wave ê incident from n̂, which solve and solve_chiral provide in every layer (Pasteur media are reciprocal, so this holds inside chiral shells). Rotating the emitter onto the z axis leaves only m = 0, ±1 multipoles, so any position and complex orientation costs one set of π_l, τ_l; total and helicity-resolved powers are summed from orthogonal multipole amplitudes. Every emitter is six basis sources (unit p and m along x, y, z), so fixed sources and orientation averages are one bilinear form in the source covariance ⟨ss†⟩, exact rather than sampled. For an emitter in the host the direct dipole field is added in closed form, so the truncation depends on the particle and not on the emitter's distance.
• Chiral shells (chiral.py). In a Pasteur medium the field splits into Beltrami waves of helicity ±1 with wavenumbers k₀(n ± κ), and an interface mixes them only through its impedance contrast. The regular solution is swept outwards as a 2 × 2 matrix ρ with the same scaling as the achiral solver: the matrix Möbius step is assembled from the index, chirality and impedance contrasts (no cancellation of the large l/x parts), and ρ is carried as element-wise logarithms so it survives falling below the double range across thick absorbing shells. The layer amplitudes follow inwards from the host through the inverse of the same step, in element-wise logarithms as well.
• Decay rates of any source (rates.py). In the emitter's layer a source (p, m, Q) drives each helicity channel s through q_s = p + i s m/Z plus (1/6)Q:∇; the direct field is the dyadic expansion of Tai, the reflected one follows from the 2 × 2 R, S of the layer evaluated scaled at the emitter, B = (1 − RS)⁻¹R(D_ψ + S D_ξ), A = S(B + D_ξ). The total rate is the free power (closed form) plus Re q*·(Aψ + Bξ); the radiated power, per helicity, comes from the outgoing amplitudes carried to the host by per-interface maps (ChiralSolution.log_in, log_out), so no ill-conditioned amplitude matrix is ever inverted; the absorption of each layer is the loss density ε″|E|² + μ″|H|² − 2κ″ Im(E*·H) integrated by Gauss–Legendre quadrature in the basis (ψ₊, ψ₋, ξ₊, ξ₋) with per-function logarithmic scales, and that of each sheet comes from the fields on its two sides. The three are independent, and balance_error checks energy conservation. Orientation averages are exact: the analytic ⟨ss†⟩ for dipoles, the 60 rotations of the icosahedral group (exact to rank 4) once a quadrupole is present.
• 2D sheets (sheets.py). The averaged sheet conditions are a unimodular 2 × 2 transfer on each pair of tangential components (E_M, H_N) and (E_N, H_M); they enter the Möbius steps of both solvers as additions to the interface matrices (the large l/x parts still cancel analytically in the bulk contrast), and the amplitude equations acquire the matching terms.
• Quadrupoles. With the emitter on the local z axis a quadrupole couples to the axial families m = 0, ±1, ±2; the gradients of M and N there are closed forms in ψ'/ψ, x and l. Their far field is synthesised from the outgoing host amplitudes (for dipoles this route agrees with the reciprocity route to 1e-14 in complex amplitude).
• Ohmic loss from boundary terms. The absorption in a shell, ∫|A j_l + B h_l|² r² dr, is a Lommel integral: it needs the fields only at the shell's boundaries, so its cost is O(l) rather than the O(l²) of quadrature (whose node count must grow with the order). Past the turning point the textbook boundary term Im(k* f_{l−1}* f_l) is a small difference of real numbers, losing ~log₁₀(2l²/Im x²) digits; the jj and hh parts are therefore carried by ratio recurrences (downward for j, upward for h) that add only positive terms there, and the jh part has no cancellation. Weakly lossy shells (|Im k| < 10⁻³|k|), where Im k² → 0, use quadrature, in chunks so memory stays bounded. An emitter 1 nm from a 1-µm silver sphere needs ~17000 orders and takes 0.8 s (quadrature ran out of memory); the default order cap is 20000.
• Incident field in closed form. In the host, near_field and energy_density sum only the scattered series, which converges with the sphere's orders at any distance, and add the plane wave exactly (its own series needs l ~ kr: with Wiscombe's truncation |E|² eight radii from a gold sphere was off by a factor of two). near_field(..., incident=False) returns the scattered field alone. A partial-wave field, near_field(..., orders=, polarisations=), is the exception: it sums the incident wave's own terms of those orders, which is what such a field is, so that disjoint selections add up to their union and all orders together to the field itself. energy_density and shell_energy take the same selection; averaged over a sphere the partial waves do not interfere, so there the selections add up in energy, not only in field.
• Exact field at the centre. At r = 0 only the l = 1 internal waves survive, so near_field evaluates the origin in closed form rather than at a small offset: E(0) = d₁ x̂ and H(0) = m c₁ ŷ for a homogeneous sphere (Bohren & Huffman), to ~1e-15 (benchmarks/MIE_CODES.md). scattnlay, which evaluates the series there, is off by up to 5e-3. scattering_amplitudes(..., orders=, polarisations=) keeps chosen partial waves of the far field: their amplitudes add up, their |S|² does not.
• Vectorised. solve takes a whole spectrum at once (array operations over wavelength × order; Python loops only over interfaces). Near fields and energy densities evaluate the special functions once per distinct kr for both polarisations.

Accuracy you can expect. Coefficients agree with 60-digit transfer matrices to ~1e-12, and the interface conditions hold to 1e-9 at every order. Where a result is less accurate than that, the problem is ill-conditioned rather than the algorithm: large lossless particles with sharp internal resonances can amplify input rounding by 10⁶ (perturbing the inputs by 1e-15 moves the exact answer by as much as the solver's error), and the absorption of a nearly lossless Rayleigh particle is known only to ~eps·|a_l|. Layers index-matched to within δ determine their coefficients to ~l·eps/δ. Chiral T-matrices agree with 80-digit 4 × 4 transfer matrices to 1e-12 up to l = 160 (the amplitudes in every layer to 6e-13), with the independent code treams to 1e-13, and a small chiral sphere reproduces the published quasi-static polarisabilities α_EE, α_HH, α_EH of Guzatov & Klimov (2012).

emission_rates equals decay_rates to ~1e-13 (total, radiative, nonradiative) for achiral spheres; for lossless chiral multilayers its total rate equals the independent reciprocity-based radiated power to 1e-11, per helicity to 1e-12; with absorbing, magnetic and chiral layers, gold and sheets, total = radiative + absorbed holds to ~1e-12 even 0.1 nm from an interface; the orientation-averaged rate of a chiral molecule beside a small chiral sphere converges to Guzatov & Klimov's quasi-static Eq. 46 as (k₀a)² for both enantiomers. On the logarithmic route (route="log"), within ~1 nm of a lossless interface the reflected part of the total rate is the small real part of large evanescent terms, so rounding limits total to ~1e-9 (1 nm) … 1e-5 (0.1 nm) relative (the route flags it through the energy balance); the normalized route, the default, keeps full precision there (see Normalized formulation below). Quadrupole couplings agree with 40-digit finite differences of the vector wave functions to 1e-14, and a pair of opposite dipoles ±p at r₀ ± δ/2, computed by the reciprocity route, converges to the quadrupole 3(pδ + δp) − 2(p·δ)I plus the current loop (ik₀/2)p × δ to 1e-12 in amplitude, near and inside chiral particles. A sheet is the O(d²)-accurate limit of the explicit film it replaces, and graphene-coated-sphere resonances sit on the quasi-static condition ε₁l + ε₂(l+1) + iσl(l+1)/(k₀R) = 0 within the (k₀R)² retardation shift.

Normalized formulation (the default route). Every per-order term at a source is an explicit geometric factor times quantities of order unity built from ψ'/ψ, ψξ and the normalized Bessel function j̄ (normalized.py): the reflection ratios ρ, σ are swept through Möbius maps whose mismatches come from the material contrasts, and near a lossless interface the small real parts of the nearly reactive high-order terms are carried multiplicatively, so total rate, radiative rate and frequency shift keep full relative precision at any distance (energy balance ≤ 1e-14 down to 0.05 nm from lossless interfaces). Interface responses - sheets (averaged transition conditions) and Feibelman d-parameters (linearized from the metal side) - add terms to the same maps through one matching form, value′ = τ(value + t_v deriv)/c_v, deriv′ = τ(deriv − t_d value)/c_d. Chiral layers enter as 2 × 2 maps in the (TM, TE) basis, with the leading mismatches written through a closed-form contrast matrix (its TE element is exactly zero between achiral nonmagnetic media, where the helicity basis would lose ε(l/x)² of the TE terms). Absorption is computed independently

• Lommel boundary terms for layers, the jump of the radial flux for interface responses - so the energy balance remains a test.

Known limitations.

• Absorption next to metal. The absorbed power (per layer and per interface response, i.e. EmissionRates.absorption, .sheet_absorption and DecayRates.nonradiative) is computed from the amplitudes of the logarithmic solver, whose logarithms grow like l ln(l/x). Within a few tenths of a nanometre of an absorbing layer, where 10³–10⁴ orders are summed, this limits the absorption to ~1e-13 relative at 1 nm, ~1e-12 at 0.1 nm and ~1e-11 at 0.05 nm (visible as the energy balance total − radiative − absorbed). Total rate, radiative rate and frequency shift come from the normalized formulation and are not affected.
• Quadrupoles deep in a subwavelength lossless core. There the outgoing solution is nearly a standing wave and the radiative part 1 + Re S of the low orders is a small difference: ~1e-14 for dipoles and 5–8e-14 for quadrupoles in a core of index contrast 3.5, worse for smaller cores and higher multipoles.
• Chiral layers. The normalized route needs every source in an achiral layer and no sheets in a structure with chiral layers; otherwise route="auto" takes the logarithmic route, with its limit next to lossless interfaces. With absorbing layers in a chiral structure the absorption is still integrated by quadrature, O(L^2): within ~0.5 nm of an absorbing layer (10^3-10^4 orders) this becomes slow (minutes to hours per position), and its accuracy is ~1e-12 at 1 nm.
• d-parameters. The mesoscopic boundary conditions are first order in d; they are applied linearized from the metal side (exact for d∥ = 0, and exactly the sphere Mie coefficients of Gonçalves et al., Nat. Commun. 11, 366 (2020)); with both d⊥ and d∥ nonzero, other second-order readings differ by d⊥d∥ l(l+1)/R².

Typical timings on one core:

Benchmarks

python benchmarks/convergence.py reruns the test cases of the multilayered-sphere literature against independent references and writes the report in benchmarks/RESULTS.md: decay rates near silver and silicon spheres down to 0.25-nm gaps (Majic & Le Ru 2020) and inside and outside gold nanoshells and a matryoshka (against high-precision transfer matrices), a dipole in a layered magnetic sphere and a Luneburg lens (Yuan, Zhu & Zhu 2023, 2024), 1000-layer spheres at x = 1000, a graded 250-µm droplet at the rainbow angle and the Cauchy profile (Wu & Wang 1991; Wu et al. 1997), absorbing spheres to x = 20000, and fields far from the particle.

The planar solver descends from Byrnes's tmm, so it is checked against two codes outside that lineage: python benchmarks/planar_codes.py (after pip install -e ".[benchmarks]") compares it with PyMoosh (scattering matrices) and pyElli (2×2 and Berreman 4×4) and writes benchmarks/PLANAR_CODES.md: R, T, r, t, absorption per layer, Ψ and Δ agree to 1e-14 in power and 6e-13° in angle, through opaque and 20-µm absorbing layers, total internal reflection, a magnetic layer and an incoherent substrate; plane-wave fields agree with PyMoosh to 2e-14, and the guided modes behind a film emitter's guided channel with PyMoosh's mode search (1e-11) and the slab dispersion relation (7e-14). python benchmarks/emitter_references.py checks film emitters: decay rates against Sommerfeld integrals written without pystratify (themselves pinned to the image dipole at a perfect mirror, to 5e-12), and near lossless guides against smuthi and the lossless limit of an absorbing guide (benchmarks/EMITTER_REFERENCES.md). python benchmarks/mnpbem_spheres.py checks spheres by a different method, the boundary elements of MNPBEM run under GNU Octave (benchmarks/mnpbem_octave.py fetches the CPC program, verifies it and builds an Octave-ready copy): cross sections and decay rates converge toward PyStratify as the mesh is refined, to 0.5–4% at the finest mesh (benchmarks/MNPBEM_SPHERES.md). python benchmarks/mie_codes.py compares spheres with two independent Mie codes, scattnlay and treams: cross sections, backscattering and amplitudes of ten spheres up to x = 1000 (nanoshells, a matryoshka, a 10-layer graded sphere) and the T-matrices of chiral multilayers agree to 2e-11 once 20 orders beyond Wiscombe's truncation are kept, and to 4e-8 at the default truncation (benchmarks/MIE_CODES.md). python benchmarks/gdm_spheres.py adds a third method, the volume Green dyadic method of pyGDM2: dielectric spheres agree within 3% in cross sections and 0.2% in decay rates 40 nm from the surface, at the volume-equivalent radius (benchmarks/GDM_SPHERES.md). Every comparison is a script in benchmarks/ that writes its own report, so it can be rerun by anyone; the tests import them. benchmarks/MATRIX.md maps every check, existing and planned, by geometry, source and quantity.

Tests

pip install -e ".[test]"
pytest

• test_planar_emission_references.py: film emitter rates (electric and magnetic, perpendicular and parallel) against the Sommerfeld reference of benchmarks/emitter_references.py, itself checked against the image dipole at a perfect mirror; sources outside lossless guides against the lossless limit and smuthi (the latter skipped unless smuthi is installed).
• test_gdm_spheres.py: pyGDM2's volume dipoles within a few percent of PyStratify for dielectric spheres (skipped unless the benchmarks extra is installed).
• test_mie_codes.py: spheres against scattnlay and treams (skipped unless the benchmarks extra is installed; the CI job benchmarks installs it).
• test_mnpbem_spheres.py: MNPBEM's boundary elements approach PyStratify's sphere cross sections and decay rates as the mesh is refined (skipped unless GNU Octave is installed).
• test_planar_codes.py: the planar solver, its fields and its guided modes against PyMoosh, pyElli and the slab dispersion relation (skipped unless the benchmarks extra is installed; the CI job benchmarks installs it).
• test_api_contract.py: every name in __all__ has exactly one place in docs/API.md (frozen, provisional, legacy or internal).
• test_focused.py: focused illumination (above): every semi-analytic form against plane-wave superpositions through the existing engines, energy balance, limits, Köhler and Parseval.
• test_normalized.py: the normalized formulation (normalized.py: ψ'/ψ, ψξ and normalized j̄ only, reflection ratios swept outwards and inwards) against the solver, extended precision and the unnormalized formulas, including emitters and interfaces at real zeros of ψ_l, the small real part Re ψξ = ψ² and strongly amplifying shells.
• test_routes.py: the normalized route of decay_rates and emission_rates against the logarithmic one for every source (dipoles, chiral p ∥ m, general p + m, all electric- and magnetic-quadrupole components, orientation averages), energy conservation 1 to 0.05 nm from lossless interfaces, the energy-balance flag of the logarithmic route, the shift of a chiral emitter as the weighted sum of its electric and magnetic shifts, and the magnetic quadrupole by duality (a dual multilayer gives it the rates of the electric one).
• test_shift.py: the per-order Green's forms of normalized_terms and the frequency shift against mpmath transfer matrices (references.layered_green_forms, layered_green_sums), the quasi-static limit −Re r/(2 Im r) of the shift near gold, the automatic truncation (which includes the reactive tail the shift needs) and the vectorized emitter side.
• test_riccati.py: special functions against mpmath in every regime above, and the Wronskian for |x| up to 8·10⁴.
• test_solver.py: coefficients against 60-digit transfer matrices (including a 1-nm film, index-matched layers and a gain shell), the interface conditions at every order up to 400, batching, 50 shells, x = 1000, Im x ≈ 700, and BHMIE at x = 20000.
• test_physics.py: textbook Mie (SciPy's spherical Bessel functions), the Bohren–Huffman reference case, the optical theorem, energy conservation of lossless multilayers, passivity, the Rayleigh limit, plane-wave limits, field continuity and convergence far from the particle, quadrature, decay rates against the mpmath Mie sums of Majic & Le Ru and, for emitters inside and outside metal shells, against high-precision transfer matrices, closed-form against quadrature Ohmic loss, magnetic-dipole rates against transfer matrices (including μ ≠ 1 and magnetically lossy shells), and energy conservation of decay rates for electric and magnetic dipoles, up to 17000 orders. Partial-wave selections of the near field, the far field and the energy add up to the whole, and an electric dipole alone has the dipole's pattern.
• test_farfield.py: the amplitude matrix against textbook Mie with independently computed angular functions, the pattern as the far-zone limit of the near field (residual falling as 1/kr), Bohren & Huffman's backscattering value, Mueller-matrix identities, directivity and dσ/dΩ integrating to 4π and σ_sca. Dipoles: closed forms in a homogeneous medium for any position and complex orientation, reciprocity against the plane-wave near field for emitters in the core, in an absorbing shell, in a magnetic shell and in the host (every phase), radiated power equal to decay_rates to 1e-12, a distant dipole reproducing the plane-wave amplitude of a chiral sphere, and helicity conservation of a (p, −ip) source beside a dual chiral particle.
• test_chiral.py: 80-digit 4 × 4 transfer matrices (thin Au shell, gain, a 1e-6 chirality contrast, l up to 160), treams, κ = 0 against the achiral solver, energy conservation of lossless chiral multilayers, reciprocity (symmetric TM–TE block), mirror symmetry (κ → −κ exchanges the helicities), helicity conservation of dual particles, and a chiral core behind 2 µm of metal, the published small-sphere polarisabilities, a chiral metal with Re(n − κ) < 0.
• test_emitters.py: emitters inside chiral particles against the achiral solver (κ → 0), the Pasteur interface conditions (tangential E and H, D_r = εE_r + iκH_r, B_r = μH_r − iκE_r) seen by emitters on either side of every interface, a (q, −iZq) source inside a dual chiral particle emitting one helicity only, the free chiral emitter's g_lum in a medium, orientation averages against the exact averages over the 24 rotations of the cube and 8 about an axis, quadrature of total and helicity-resolved power, and the mirror symmetry of valley excitons beside spheres (g_lum = 0 alone, ±g on either side or valley).
• test_rates.py: decay_rates (κ → 0), the reciprocity far field for lossless chiral multilayers (total = radiated, per helicity), energy balance with absorbing chiral, magnetic and gold layers down to 0.1 nm from interfaces, mirror symmetry (κ, m) → (−κ, −m), exact orientation averages, batches of positions, the free power in a chiral medium, normalisations and quantum yield, and Guzatov & Klimov's Eq. 46 for dielectric, lossy, metallic and double-negative chiral spheres.
• test_sheets.py: explicit films (O(d²) with the out-of-plane term, O(d) without it), sheets between chiral layers, the chiral and achiral solvers against each other to l = 400, graphene-coated-sphere plasmons against the quasi-static condition, graphene conductivity limits, plane-wave absorption = Poynting-flux jump = sheet-loss formula = extinction − scattering, emitter energy balance and far fields with sheets.
• test_quadrupoles.py: the coupling against 40-digit finite differences, Jackson's free power, total = radiated and energy balance with quadrupoles, the dipole-pair limit in amplitude, the two far-field routes against each other, pattern quadrature, the icosahedral design, and current-loop moments.
• test_pairs.py: the two-point Green's dyadic (green_dyadic): the multipole series of the free field against the closed-form free dyadic, reciprocity, the coinciding limit (decay rates and shift), the radial two-point forms in one shell and across a gold shell (transmission chain) against mpmath order by order, and the interface conditions approached from either side.
• test_feibelman.py: Feibelman d-parameters against the mesoscopic Mie theory of Gonçalves et al. (2020), Eqs. (4a) and (13a); the per-order forms of a nanoshell with d-parameters on both surfaces against mpmath matrices built from the boundary conditions; energy balance including the surface absorption; the TE response as the sheet σ = ik₀d∥(ε_other − ε_metal).
• test_chiral_normalized.py: chiral layers in the normalized formulation: the achiral limit is the scalar sweep, the 2 × 2 reflection matrices against 80-digit 4 × 4 transfer matrices, energy balance, agreement with the logarithmic route and mirror symmetry for chiral emitters, and emitters inside a chiral layer routed to the logarithmic formulation.
• test_ensemble.py: shell averages (free space, graded nodes, the d_min⁻² cutoff next to gold, the thin-shell limit) and the spectral density J(ω).
• test_casimir.py: the Casimir–Polder sums at imaginary frequency against mpmath transfer matrices, the non-retarded limit −C₃/d³ next to a large Drude sphere, and the convergence of the frequency quadrature.

pystratify.references holds the independent references (mpmath Mie decay rates and extinction, layered-sphere decay rates from transfer matrices, BHMIE); nothing in it shares code with the package.

Citing

PyStratify implements the method of STRATIFY; please cite the paper,

I. L. Rasskazov, P. S. Carney and A. Moroz, "STRATIFY: a comprehensive and versatile MATLAB code for a multilayered sphere," OSA Continuum 3, 2290 (2020), doi:10.1364/OSAC.399979,

and the software itself: doi:10.5281/zenodo.23287561 names every version and resolves to the latest; each release also has its own DOI on that Zenodo page, which is the one to cite for the exact code you ran (CITATION.cff). Defects of the MATLAB code that PyStratify does not share are listed in AUDIT.md.

Licence

GPL-3.0-or-later, as STRATIFY (this is a derivative work).