treaaams v0.4.6c6

A derivative of treams with targeted source-level patches to enable complex frequency evaluation.

What treaaams changes in treams

treams uses the time convention $\mathrm{Re}(\boldsymbol E(\boldsymbol r, k_0),\mathrm e^{-\mathrm i c k_0 t})$ (treams paper, Sec. 2.1), so resonances (poles of T- and S-matrices) lie at $\mathrm{Im},k_0 < 0$. To reach them, every quantity has to be the analytic continuation of its value at real $k_0$. treams instead picks branches pointwise, e.g. "the root with $\mathrm{Im} \geq 0$", which is correct on the real axis but jumps when $k_0$ crosses it. treaaams makes the changes below. At real $k_0$ the physical results equal those of treams. The upstream test suite passes except for two tests that evaluate $\Gamma(n, z)$ directly on its moved branch cut (see Tests).

1. Complex $k_0$ is accepted everywhere. TMatrix, TMatrixC, SMatrix and the k0 annotation accept complex values. The basis CylindricalWaveBasis accepts complex $k_z$. The Mie coefficients of spheres take complex size parameters $x = k_0 r$, and those of cylinders take complex $k_0$ and $k_z$ (patches 020–024, 050).

Branch cut convention. Square roots, powers and $\Gamma(n, \cdot)$ are multivalued. treams evaluates them on their principal branches, with the cut along the negative real axis, and in places adjusts the sign afterwards (see 2.) (DLMF 4.2.5, 4.2.28). treaaams places the cut of each of them along a configurable ray $\arg\zeta = \varphi$, with $\varphi \in [-\pi, \pi]$. Of the two branches bounded by that ray, it takes the one containing the positive real axis: $\arg\zeta$ is taken in $(\varphi - 2\pi, \varphi]$ if $\varphi \geq 0$, and in $(\varphi, \varphi + 2\pi]$ if $\varphi < 0$. Thus $\varphi = \pm\pi$ is the principal branch, and all functions agree with treams for real positive arguments. The settings are listed under Configuration.

2. Square roots of the normal wave vector component are analytic. Plane waves use $k_z = \sqrt{\zeta}$ with $\zeta = k^2 - k_x^2 - k_y^2$, and cylindrical waves use $k_\rho = \sqrt{k^2 - k_z^2}$. treams takes the root with $\mathrm{Im} \geq 0$ (see treams paper, text after Eq. (E.1), and lattice paper, text after Eq. (19)). treaaams uses the cut $\varphi =$ BRANCH_CUT_SQRT_KZ $= -\pi/2$ instead (misc.angled_sqrt, patches 023, 030–032, 041):

$$ \sqrt{\zeta} = \begin{cases} -\sqrt{\zeta}{\mathrm{p}} & \arg \zeta \leq -\pi/2 \ \phantom{-}\sqrt{\zeta}{\mathrm{p}} & \text{otherwise} \end{cases} \qquad\Rightarrow\qquad \arg\sqrt{\zeta} \in \left(-\tfrac{\pi}{4}, \tfrac{3\pi}{4}\right], $$

where $\sqrt{\zeta}{\mathrm p}$ is the principal root (DLMF 4.2.28). For real $\zeta$ this gives the same root as treams. Unlike treams, though, propagating orders ($\zeta$ close to the positive real axis) and evanescent orders ($\zeta$ close to the negative real axis) are continued smoothly into $\mathrm{Im},k_0 < 0$. There, the propagating orders grow away from the scatterer, as quasi-normal modes do. All later sign adjustments of $k_z$ and $k\rho$ (in sw, cw, the field evaluations and cw.translate_periodic) are removed.

3. The incomplete gamma function of the Ewald lattice sums is continued.

Where it appears. The reciprocal-space parts of the lattice sums contain the upper incomplete gamma function $\Gamma(n, z) = \int_z^\infty t^{n-1}\mathrm e^{-t},\mathrm dt$ (DLMF 8.2.2) at $z = -\gamma^2/(2\eta^2)$, with $\gamma = \sqrt{1 - \beta^2}$, $\beta = |\boldsymbol k_\parallel + \boldsymbol G|/k$ and the Ewald parameter $\eta$:

• $\Gamma(1, z)$: lattice paper, Eq. (20)
• $\Gamma(\tfrac12 - n, z)$: lattice paper, Eqs. (25) and (35)
• $\Gamma(-n, z)$: lattice paper, Eq. (30)
• origin terms $\Gamma(0, -\tfrac{1}{2\eta^2})$ and $\Gamma(-\tfrac12, -\tfrac{1}{2\eta^2})$: lattice paper, Eq. (14)

The problem. For a propagating order at real $k_0$, $\beta < 1$, so $z$ is negative real. That is exactly the cut of the principal $\Gamma(n, z)$, and the lattice paper prescribes the value from below, $z - \mathrm i 0$ (text after Eq. (14)). treams does the same. Once $k_0$ is complex, $z$ moves off the axis. Below it, the principal value continues the value from below, but above it the principal value belongs to the other side of the cut and jumps.

The fix. treaaams moves the cut out of the way, to the positive imaginary axis ($\varphi =$ BRANCH_CUT_INCGAMMA $= \pi/2$, patches 011–014). Points between the new and the old cut ($\pi/2 < \arg z \leq \pi$, second quadrant) are no longer evaluated on the principal sheet. They are evaluated on the sheet reached by continuing from below across the negative real axis. Following DLMF, that is written $\Gamma(n, z,\mathrm e^{-2\pi\mathrm i})$: the same point $z$, reached by one clockwise turn around the origin. Everywhere else the principal value is kept. The function with the original cut is still available as special.unshifted_incgamma.

The value on the neighbouring sheet depends on the order $n$; with $m = -1$ in this case:

order $n$	$\Gamma(n, z,\mathrm e^{2\pi\mathrm i m})$	why
$1, 2, 3, \dots$	$\Gamma(n, z)$	the integrand $t^{n-1}\mathrm e^{-t}$ has no branch point, so there is only one sheet
$\pm\tfrac12, \pm\tfrac32, \dots$	$2,\Gamma(n) - \Gamma(n, z)$	DLMF 8.2.9 with $\mathrm e^{2\pi\mathrm i m n} = -1$
$0, -1, -2, \dots$	$\Gamma(n, z) - 2\pi\mathrm i, m, \dfrac{(-1)^{n}}{(-n)!}$	derived below

For the half-integer orders, the general continuation formula DLMF 8.2.9,

$$ \Gamma(n, z,\mathrm e^{2\pi\mathrm i m}) = \mathrm e^{2\pi\mathrm i m n},\Gamma(n, z) + \left(1 - \mathrm e^{2\pi\mathrm i m n}\right)\Gamma(n), $$

applies directly. For $n = -N$ with $N = 0, 1, 2, \dots$ it can't be used, since $\Gamma(-N)$ is infinite while $1 - \mathrm e^{2\pi\mathrm i m n} = 0$. Instead, the last row follows in three steps:

1. Split off the multivalued part. DLMF 8.4.15 writes $$\Gamma(-N, z) = \frac{(-1)^N}{N!}\left(E_1(z) - \mathrm e^{-z}\sum_{k=0}^{N-1}\frac{(-1)^k k!}{z^{k+1}}\right).$$ The sum contains only integer powers of $z$, which are the same on every sheet. All of the multivaluedness sits in the exponential integral $E_1(z) = \Gamma(0, z)$ (DLMF 8.4.4).
2. Continue $E_1$. $E_1$ has a logarithmic branch point at $0$, like $-\ln z$. Each turn around the origin therefore shifts it by a constant: $E_1(z,\mathrm e^{2\pi\mathrm i m}) = E_1(z) - 2\pi\mathrm i m$ (DLMF 6.4.2).
3. Insert. Only the $E_1$ term changes, so $$\Gamma(-N, z,\mathrm e^{2\pi\mathrm i m}) = \Gamma(-N, z) - 2\pi\mathrm i, m,\frac{(-1)^N}{N!}.$$

As a check, for $N = 0$ and $m = -1$ the shift is $+2\pi\mathrm i$. That is exactly the jump of $E_1$ across the negative real axis, from $E_1(-x + \mathrm i 0)$ to $E_1(-x - \mathrm i 0)$ (DLMF 6.4.3), which the moved cut removes.

The accompanying power. The sums use the reduced function $\Gamma(n, z)/(-z)^n$ (_esum._redincgamma, patch 015). Its power $(-z)^n = \exp(n \operatorname{Ln}(-z))$ (DLMF 4.2.26) uses the cut $\arg(-z) =$ BRANCH_CUT_CPOW $= -\pi/2$, which is the same ray as $\arg z = \pi/2$. The regularization near $z = 0$ can be configured through SINGULARITY_REDINCGAMMA and SINGULARITY_THRESH_REDINCGAMMA.

4. The direction of plane waves is given explicitly in the lattice expansions. Periodic arrays of spherical waves on a 2D lattice expand into plane waves with the factor $\left(A k^2 \sqrt{1 - (\boldsymbol k_\parallel + \boldsymbol G)^2/k^2}\right)^{-1} = (A k k_z)^{-1}$ (treams paper, Eq. (35)). Arrays of cylindrical waves expand with $\propto (a k_y)^{-1}$ (Eq. (36)). The sign of the normal component follows the side of the lattice, $k_z = \pm k\gamma$ for $z \gtrless 0$ (text after Eq. (E.1)). treams applies that sign through a sign adjustment to $\mathrm{Im},k_z \geq 0$, which is not analytic. In treaaams, sw.periodic_to_pw and cw.periodic_to_pw therefore use $d,k_z$ (or $d,k_y$) with an explicit direction $d = \pm 1$ (patches 030, 040, 042, 043). The operators pass $d$ from the basis ($+1$ for up and $-1$ for down). For bases given by unit vectors, $d$ is the sign of the unit vector's normal component. The keyword direction is optional. Without it, $d$ is the treams sign of $k_z$ (or $k_y$), so results match treams, but they are then not analytic.

5. Mie coefficients. The refractive index $n = \sqrt{\epsilon\mu}$ and the relative impedance $Z = \sqrt{\mu/\epsilon}$ (treams paper, text after Eq. (4)) take their roots with configurable cuts, BRANCH_CUT_SQRT_MIE_N and BRANCH_CUT_SQRT_MIE_Z (default $-\pi$, i.e. the principal root as in treams). The same root gives $k_\rho$ in the cylinder coefficients (patch 050).

Other changes. The package is renamed to treaaams, so it can be installed next to treams. BLAS is limited to one thread on import (010). The warning for mismatched positions and T-matrices now shows both counts (090).

Configuration

Every cut follows the branch cut convention and can be read and set in treaaams.config with get_<NAME>() and set_<NAME>(value):

name	default	cut along
BRANCH_CUT_SQRT_KZ	$-\pi/2$	$\arg(k^2 - k_\parallel^2) =$ value
BRANCH_CUT_INCGAMMA	$\pi/2$	$\arg z =$ value
BRANCH_CUT_CPOW	$-\pi/2$	$\arg(-z) =$ value
BRANCH_CUT_SQRT_MIE_N, BRANCH_CUT_SQRT_MIE_Z	$-\pi$	$\arg(\epsilon\mu) =$ value, $\arg(\mu/\epsilon) =$ value

References:

• [treams paper]: D. Beutel, I. Fernandez-Corbaton, C. Rockstuhl, treams – A T-matrix scattering code for nanophotonic computations, arXiv:2309.03182 (2023). Equation numbers refer to v1.
• [lattice paper]: D. Beutel, I. Fernandez-Corbaton, C. Rockstuhl, Unified lattice sums accommodating multiple sublattices for solutions of the Helmholtz equation in two and three dimensions, Phys. Rev. A 107, 013508 (2023), doi:10.1103/PhysRevA.107.013508, arXiv:2207.06930. Equation numbers refer to the arXiv version.
• DLMF: NIST Digital Library of Mathematical Functions, §§ 4.2, 6.4, 8.2, 8.4.

How it works

GitHub (treams vX.Y.Z)
        │  tarball download
        ▼
_build/treams_src/          ← clean upstream source
        │  patch -p1 (each file in patches/ in order)
        ▼
_build/treams_src/          ← patched source
        │  build (and publish) wheel
        ▼
install treaaams via pip/uv

Prerequisites

• Python ≥ 3.10
• uv (recommended) or pip
• A C compiler + patch utility (apt install patch / brew install gpatch)
• Cython and NumPy (installed automatically by the build backend)

Quick start: editable install

# 1. Create a virtual environment (uv)
uv sync

# 2. Fetch upstream, apply patches, and install in editable mode
make editable

# 3. Verify
python -c "import treaaams; print(treaaams.__version__)"

Tests

make test           # tests in tests/
make test-upstream  # upstream treams unit tests run against treaaams

Upstream tests that intentionally differ on treaaams are marked as strict xfail in tests/upstream_xfail.py, together with the reason. CI (.github/workflows/tests.yml) runs both on every push to main and on every pull request for Python 3.10–3.13, and publishing a tag only proceeds when they pass.

Adding or modifying a patch

Workflow

# 1. Make sure you have a clean source tree
python scripts/build_patched.py --no-fetch   # skips download, resets patches

# 2. Edit the file you want to change
$EDITOR _build/treaaams_src/src/treams/_special/sph.pyx

# 3. Capture the diff as a numbered patch file
python scripts/make_patch.py src/treaaams/_special/sph.pyx \
    --name 010-custom-sph

# 4. Inspect it
cat patches/010-custom-sph.patch

# 5. Dry-run all patches to verify
make patch-dry-run

# 6. Rebuild
make editable

Patch naming convention

Patches are applied lexicographically, so number them explicitly:

File	Applied
000-fix-typo.patch	first
010-override-ufunc.patch	second
020-custom-lattice.patch	third

Tracking a new upstream version

Run for example

uv run tbump 0.4.7c1

To apply the patches to v0.4.7 of treams

If a patch no longer applies cleanly, patch will create a .rej file with the rejected hunks. Edit the patch accordingly and re-run.