# Computing the T-matrix of arbitrarily shaped scatterers: addatmatrix wrapper for the ADDA code

Maxim A. Yurkin, Krzysztof M. Czajkowski, Lorenzo Pattelli

► To cite this version:

Maxim A. Yurkin, Krzysztof M. Czajkowski, Lorenzo Pattelli. Computing the T-matrix of arbitrarily shaped scatterers: addatmatrix wrapper for the ADDA code. Bremen Workshop on Light Scattering 2025, University of Bremen, Mar 2025, Bremen, Germany. pp.17-22. hal-05060826

HAL Id: hal-05060826

https://hal.science/hal-05060826v1

Submitted on 8 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Computing the T-matrix of arbitrarily
shaped scatterers: **addatmatrix**
wrapper for the ADDA code

Maxim A. Yurkin¹, Krzysztof M. Czajkowski², and
Lorenzo Pattelli³

¹Université Rouen Normandie, INSA Rouen Normandie, CNRS,
CORIA UMR 6614, Rouen, France, yurkin@gmail.com

²Faculty of Physics, University of Warsaw, Warsaw, Poland,
Krzysztof.Czajkowski@fuw.edu.pl

³Istituto Nazionale di Ricerca Metrologica, Turin, Italy,
l.pattelli@inrim.it

The discrete dipole approximation [1] is a general method to simulate light scattering by particles with arbitrary shape and composition. It is highly optimized for a single wave incident on a particle; thus, multiple independent simulations are required if this incident wave or particle orientation are varied. By contrast, the T-matrix of the particle encodes the response of a particle to any incident wave and allows analytical orientation averaging [2]. Historically, the T-matrix has been mostly computed using surface-discretization methods, such as extended boundary-condition or null-field methods, however, any simulation method, including the DDA, can be used to compute the T-matrix. A straightforward approach is to illuminate the particle with vector spherical wave functions (VSWFs) one at a time and to expand the scattered field in VSWFs [3]. While there exist more elaborate options based on the internal structure of the DDA equations [4], they do not fully employ the block-Toeplitz structure of the DDA interaction matrix leading to poorer scaling of computational time with the number of discretization voxels.

The only drawback of the VSWF illumination is that such
incident field is not included as a built-in option in any popular
open-source DDA codes (although some allow reading of inci-
dent wave from pre-generated file). A general solution for this is-
sue is to compute the T-matrix in plane-wave basis (both for in-
cident and scattered waves) and then transform it to the VSWF
basis (or vector spherical harmonics, VSHs, when the field is
considered only on a far-field sphere) [5]. We have implemented
this approach in a Python wrapper for the ADDA code. The
latter is a robust open-source implementation of the DDA [6],
optimized for high performance, supporting distributed memory
clusters, multi-core processors, and GPU acceleration.

The first goal of this work is to make T-matrix calculations
as easy to use for the community as possible. The addatmatrix
is an open-source Python code [7], which can handle installation
of ADDA itself on any platform with the GCC compiler (was
tested on both Windows and Linux environments, specifically
under Python 3.11). For those familiar with ADDA command-
line interface (CLI), the easiest is to use the corresponding CLI
of addatmatrix. For example, the following calculates the T-
matrix of a sphere with refractive index of 2, diameter of 200
nm at a wavelength of 500 nm:

./addatmatrix_cli.py tmatrix --adda_string "-shape sphere -m 2
  ↪ → 0 -size 200 -lambda 500"

Alternatively, one can use **addatmatrix** as a Python package
for integration into other Python-based projects. Moreover, the
resulting T-matrix is stored accordingly to a recently proposed
standardized data format [8], for which there already exist a
number of postprocessing tools. A mini-tutorial on using the
wrapper will be presented at the conference; further examples
are presented at its repository.

The second goal is to minimize the number of ADDA runs,
since that is the bottleneck of these simulations. Note that the
scaling with particle size parameter *x* ≥ 10 is unfavorable. The
time of one DDA simulation scales at least as *O*(*x*⁴), assuming
that number of iterations to solve the linear system is *O*(*x*) – it
can be much worse for large real refractive index. Moreover, the
simulation of scattered fields at *O*(*x*²) directions (corresponding
to the number of VSWFs up to the degree $L \sim x$) is $O(x^5)$, but the prefactor is much smaller so it is not expected to be a significant factor for moderate values of $x$. The number of DDA runs (incident waves) is also $O(x^2)$, resulting in the total computational time of $O(x^6)$ or greater. This is, in principle, comparable to the standard T-matrix computations based on inversion of matrices of size $O(x^2)$.

Leaving apart the optimization of the DDA itself, which is also a field of active research, we took a closer look at the plane-wave to VSWFs transformation. Mathematically, this boils down to spherical harmonic transform (SHT) [9] with some complications from generalizations to vector harmonics. The most straightforward approach to the SHT is based on least-square solution (fitting the coefficients to reproduce the data), and it is also most widely used in light-scattering community [5]. The drawbacks is that the procedure involves inversion of matrices with size $O(L^2)$, which condition number rapidly grows with $L$. Hence, while theoretically the number of data samples can be equal to the number of scalar or vector spherical harmonics, $(L + 1)^2$ or $2L(L + 2)$, respectively, it often needs to be taken significantly larger to improve the accuracy.

One of the possible improvements is the use of spherical quadratures, which set the position of the points and weights to exactly integrate the spherical harmonics up to a given degree. Due to orthogonality of VSHs the expansion coefficients of any function on a sphere can be obtained by taking its inner product with the corresponding VSHs. The integrand is a dot product of two VSHs (with indices $l, m$ and $l', m'$, the second one is conjugated), which can be shown to be a sum of products of associated Legendre polynomials (of argument $\cos\theta$) of degrees $l, l'$ (or smaller) times $\exp[i(m - m')\phi]$. This can be shown (e.g., by the Gaunt's formula) to be always orthogonal to any spherical harmonics of degree larger than $2L$. The latter is, thus, the required degree of an appropriate quadrature. For example, it can be a tensor product of Gauss-Legendre quadrature with $L + 1$ samples of $\cos\theta$ and trapezoidal rule with $2L + 1$ samples over $\phi$ (so-called Gauss-product quadrature). This approach has
been used in [10] albeit with slightly larger (redundant) number
of samples. The main advantage is that no matrix needs to be
inverted ensuring good numerical stability. However, the num-
ber of DDA simulations is almost twice larger than the number
of VSHs, since two incident polarization need to be computed
for each incident direction. While this number can be reduced
by the use of more cumbersome Lebedev quadratures [11], it will
still be at least 4/3 times larger than the optimal.

Another approach is an optimal-dimensionality sampling [9]
that can be considered a compromise between the previous two.
The grid is structured into iso-latitude rings (as for the Gauss-
product quadrature), but the number of $\phi$ samples in each ring
is different. The rings are placed uniformly in $\theta$, but their order
(with increasing number of $\phi$ samples) is optimized to minimize
the condition number of the matrix to be inverted. This condi-
tion number is then within 100 for $L \le 250$. Importantly the
total number of samples to perform scalar SHT is $(L+1)^2$ –
the minimum possible number. The main drawback is that this
algorithm needs to be extended to vector SHT.

The current version of addatmatrix, as discussed in [8],
is based on the least-square SHT, but we are working to im-
plement Gauss-product quadrature and optimal-dimensionality
sampling. The results will be presented at the conference.

To conclude, addatmatrix is an easy-to-use Python package
enabling computation of T-matrix of arbitrarily shaped inho-
mogeneous particles with the DDA. While it does not suffer
from numerical instabilities of traditional methods to compute
T-matrix, e.g. for highly prolate or oblate shapes, its compu-
tational requirements steeply increase with particle size. Thus,
optimizations aimed at minimizing the number of independent
DDA simulations are on the way. Another future extension
is the computation of T-matrix of a particle on top or near
a semi-infinite substrate, using the corresponding functionality
of ADDA.
References

[1] M.A. Yurkin, *Discrete dipole approximation*, in Light, Plasmonics and Particles, M.P. Menguc and M. Francoeur, Eds., Elsevier, Amsterdam (2023), pp. 167–198.

[2] M.I. Mishchenko, Electromagnetic Scattering by Particles and Particle Groups: An Introduction, Cambridge University Press, Cambridge, UK (2014).

[3] V.L.Y. Loke, T.A. Nieminen, N.R. Heckenberg, and H. Rubinsztein-Dunlop, *T-matrix calculation via discrete dipole approximation, point matching and exploiting symmetry*, J. Quant. Spectrosc. Radiat. Transfer **110** (2009), 1460–1471.

[4] D.W. Mackowski, *Discrete dipole moment method for calculation of the T matrix for nonspherical particles*, J. Opt. Soc. Am. A **19** (2002), 881–893.

[5] M. Fruhnert, I. Fernandez-Corbaton, V. Yannopapas, and C. Rockstuhl, *Computing the T-matrix of a scattering object with multiple plane wave illuminations*, Beilstein J. Nanotechnol. **8** (2017), 614–626.

[6] M.A. Yurkin and A.G. Hoekstra, *The discrete-dipole-approximation code ADDA: Capabilities and known limitations*, J. Quant. Spectrosc. Radiat. Transfer **112** (2011), 2234–2247.

[7] K. Czajkowski et al., *addatmatrix repository*, https://gitlab.com/k.czajkowski/addatmatrix/ (2025).

[8] N. Asadova et al., *T-matrix representation of optical scattering response: Suggestion for a data format*, J. Quant. Spectrosc. Radiat. Transfer **333** (2025), 109310.

[9] Z. Khalid, R.A. Kennedy, and J.D. McEwen, *An optimal-dimensionality sampling scheme on the sphere with fast spherical harmonic transforms*, IEEE Trans. Sign. Proc. **62** (2014), 4597–4610.
[10] M. Ganesh and S.C. Hawkins, *T-matrix computations for light scattering by penetrable particles with large aspect ratios*, J. Quant. Spectrosc. Radiat. Transfer **334** (2025), 109346.

[11] J. Burkardt, *Quadrature rules for the sphere*, https://people.sc.fsu.edu/jburkardt/datasets/sphere_lebedev_rule/sphere_lebedev_rule.html (2010).

[12] M.A. Yurkin and M. Huntemann, *Rigorous and fast discrete dipole approximation for particles near a plane interface*, J. Phys. Chem. C **119** (2015), 29088–29094.