# First results of using weighted discretization in the ADDA code

Paul Bouillon, Maxim A. Yurkin

► To cite this version:

Paul Bouillon, Maxim A. Yurkin. First results of using weighted discretization in the ADDA code. Bremen Workshop on Light Scattering 2025, University of Bremen, Mar 2025, Bremen, Germany. pp.23-27. hal-05060835

HAL Id: hal-05060835

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

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
# First results of using weighted discretization in the ADDA code

P. Bouillon¹ and M. A. Yurkin¹

¹ Université Rouen Normandie, INSA Rouen Normandie, CNRS, CORIA UMR 6614, Rouen, 76000, France, paul.bouillon@coria.fr

In order to describe scattering and absorption of electromagnetic waves by an arbitrary shaped particle [1], the discrete dipole approximation (DDA) is widely used. It consists in a discretization of the particle volume using a set of small voxels (dipoles), Figure 1a. Assuming the electric fields to be constant inside each voxel, the corresponding value $\mathbf{E}_i$ for a given voxel is determined from a linear system of equations, which takes into account the incident field $\mathbf{E}_i^{\text{inc}}$ and interactions $\overline{\mathbf{G}}_{ij}$ with other voxels of the set as

$$ \mathbf{E}_i = \mathbf{E}_i^{\text{inc}} + \sum_{j \neq i} \overline{\mathbf{G}}_{ij} V_j \chi_j \mathbf{E}_j + (\overline{\mathbf{M}}_i - \overline{\mathbf{L}}_i) \chi_i \mathbf{E}_i, \quad (1) $$

where $\chi_i$ is the susceptibility associated with the relative (to the host medium) refractive index $m$ of the scatterer while $\overline{\mathbf{L}}_i$ and $\overline{\mathbf{M}}_i$ are, respectively, the self-term dyadic and a correction associated to the finiteness of the voxel volume $V_i$. The traditional discretization scheme (TD) uses identical cubical sub-volumes as voxels. However, due to the mismatch of some voxels with the boundary of the scatterer (i.e., staircasing effect), shape errors are expected for most of the particles [2].

Thus, a weighted discretization (WD) is proposed to reduce these errors without changing the regular rectangular grid. First, the simplest option is to consider the effective medium approximation [3] to modify the susceptibility of these boundary voxels. Second, Piller [4] proposed an alternative based on the splitting of these voxels into two different parts (Figure 1b) by a plane interface: a primary part (p), which contains the voxel center $\mathbf{r}_i$, and a secondary one (s). The electric fields in both
domains are linearly related $E_i^s = \overline{T}_i E_i^p$ through the boundary-condition tensor $\overline{T}_i$. Therefore, an effective susceptibility $\overline{\chi}_i^e$ (2) is determined from a weighted average of the susceptibility, based on the volume fraction $f_i \triangleq V_i^p/V_i$ of the voxel:

$$ \overline{\chi}_i^e \triangleq f_i \chi_i^p \overline{\mathbf{I}} + (1 - f_i) \chi_i^s \overline{\mathbf{T}}_i. \qquad (2) $$

<figure><img src="image_2.png" /><figcaption>(a) Discretization of an arbitrary (b) Weighted discretization of a shaped. boundary dipole.</figcaption></figure>

<caption>Figure 1: Weighted discretization applied to the particle.</caption>

Although the self-term dyadic $\overline{\mathbf{L}}_i$ is determined as the sum of both parts as $\overline{\mathbf{L}}_i = \overline{\mathbf{L}}_i^p + \overline{\mathbf{L}}_i^s$, some cumbersome computational geometry is used [5] to evaluate it for each part. In the DDA framework, the general expression for the self-term dyadic is given by the following integral form [1]

$$ \overline{\mathbf{L}}(\partial V_i, \mathbf{r}_i) \triangleq \oint_{\partial V_i} \frac{\hat{\mathbf{n}}(\mathbf{r}') \otimes (\mathbf{r}' - \mathbf{r}_i)}{|\mathbf{r}_i - \mathbf{r}'|^3} d^2 \mathbf{r}' \qquad (3) $$

Since the boundary $\partial V_i$ of any polyhedron consists of several facets, one only needs to compute the contribution $\overline{\mathbf{L}}(S_i^i, \mathbf{r}_i)$ of the self-term dyadic over a single facet $S_i^i$ with vertices $\{\mathbf{p}_k\}$ and solid angle $\Omega$ (as viewed from $\mathbf{r}_i$). This contribution is deduced after some algebra, using $\hat{\mathbf{q}}_k \triangleq (\mathbf{p}_{k+1} - \mathbf{p}_k)/|\mathbf{p}_{k+1} - \mathbf{p}_k|$, to be [5]

$$ \overline{\mathbf{L}}(S_i^i, \mathbf{r}_i) = \Omega \hat{\mathbf{n}}_i \otimes \hat{\mathbf{n}}_i + \sum_{k=0}^{n-1} \ln \left( \frac{\mathbf{p}_{k+1} \cdot \hat{\mathbf{q}}_k + |\mathbf{p}_{k+1}|}{\mathbf{p}_k \cdot \hat{\mathbf{q}}_k + |\mathbf{p}_k|} \right) \hat{\mathbf{q}}_k \times \hat{\mathbf{n}}_i \otimes \hat{\mathbf{n}}_i. \qquad (4) $$
<figure><img src="image_3.png" /><figcaption>(a) Boundary voxel cut by a plane (b) Facet of the boundary voxel interface.

with vertices **p** and solid angle Ω.</figcaption></figure>

Figure 2: Computational geometry considerations for $\overline{\mathbf{L}}_i$, reproduced from [5].

Along with the above modifications, the term $(\overline{\mathbf{M}}_i - \overline{\mathbf{L}}_i) \chi_i \mathbf{E}_i$
in (1) becomes $\left[ (\overline{\mathbf{M}}_i^\mathrm{p} - \overline{\mathbf{L}}_i^\mathrm{p}) \chi_i^\mathrm{p} + (\overline{\mathbf{M}}_i^\mathrm{s} - \overline{\mathbf{L}}_i^\mathrm{s}) \chi_i^\mathrm{s} \overline{\mathbf{T}}_i \right] \mathbf{E}_i$. $\overline{\mathbf{M}}_i^{\mathrm{p,s}}$
are hard to evaluate, but they decrease with voxel size $d$ unless
the particle boundary comes very close to its center. Thus, they
can be either neglected or replaced by the traditional expres-
sion $\overline{\mathbf{M}}_i$ without significant loss of accuracy (since the number
of boundary voxels is relatively small).

In order to show how the weighted discretization improves
the simulation accuracy, errors with respect to the Mie theory
for both traditional and weighted discretization schemes are in-
vestigated for spheres with radius $R = \lambda/2$ and $\lambda$ (where $\lambda$ is
the incident wavelength). For that, the light scattering code
ADDA [6] is used, for which a prototype of the weighted dis-
cretization has been implemented on a separate branch [5]. This
prototype is compared to the main version of ADDA employing
traditional discretization. Considering simple spheres of refrac-
tive index $m = 1.33$ and $m = 2+i$ in vacuum, the relative errors
on extinction and scattering efficiencies were determined (using
the Lorenz–Mie theory as a reference [7]) as a function of grid
size along the sphere diameter (Figures 3 & 4).
<figure><img src="image_4.png" /><figcaption>Figure 3: Error of extinction efficiencies for non-absorbing spheres versus number of voxels per diameter in log-log scale.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Error of extinction efficiencies for absorbing spheres versus number of voxels per diameter in log-log scale.</figcaption></figure>

In some cases, WD greatly improves the accuracy of the DDA, but this improvement is not universal. That is partially expected. First, the WD addresses only the shape errors, this will not make a significant difference, when discretization errors dominate. Second, TD is more sophisticated than it seems, because ADDA applies so-called volume correction, which stretches the whole scale of the voxel set to exactly match the volume of the original particle. This is known to be especially efficient for small particles, where the total size (or volume) of the particle comes at some power as a factor into any scattering quantity.
We will provide data at the conference for a wider range of particles (*m*, *R*, and shape) and DDA variants. The latter will include the TD without volume correction and combination of various discretization regimes with the integration of Green's tensor (IGT). It has been hypothesized that coupling IGT and WD should lead to a quadratic convergence (error vs. voxel size) [2].

## References

[1] M. A. Yurkin, *Discrete dipole approximation*, in Light, Plasmonics and Particles (2023), M. P. Menguc et M. Francœur, Éd., Amsterdam: Elsevier, pp. 167–198.

[2] M. A. Yurkin, V. P. Maltsev, et A. G. Hoekstra, *Convergence of the discrete dipole approximation. I. Theoretical analysis*, J. Opt. Soc. Am. A **23** (2006), 2578–2591.

[3] Y. Zhu, C. Liu, et M. A. Yurkin, *Reducing shape errors in the discrete dipole approximation using effective media*, Opt. Express **31** (2023), 43401–43415.

[4] N. B. Piller, *Influence of the edge meshes on the accuracy of the coupled-dipole approximation*, Opt. Lett. **22** (1997), 1674–1676.

[5] A. A. Manina, *Use of weighted discretization in the discrete dipole approximation [in Russian]*, BSc thesis (2022), Novosibirsk State University, Novosibirsk, Russia.

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

[7] C. F. Bohren et D. R. Huffman, *Absorption and Scattering of Light by Small Particles* (1983). New York: Wiley.