# IFDDA, an easy-to-use code for simulating the field scattered by 3D inhomogeneous objects in a stratified medium: tutorial

Patrick C Chaumet, Daniel Sentenac, Guillaume Maire, Md Rasedujjaman,
Ting Zhang, Anne Sentenac

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

## ► To cite this version:

Patrick C Chaumet, Daniel Sentenac, Guillaume Maire, Md Rasedujjaman, Ting Zhang, et al. IFDDA, an easy-to-use code for simulating the field scattered by 3D inhomogeneous objects in a stratified medium: tutorial. Journal of the Optical Society of America. A Optics, Image Science, and Vision, 2021, 38 (12), pp.1841. 10.1364/josaa.432685. hal-03462189

HAL Id: hal-03462189

https://hal.science/hal-03462189

Submitted on 1 Dec 2021

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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# IFDDA, an easy-to-use code for simulating the field scattered by 3D inhomogeneous objects in a stratified medium: tutorial

PATRICK C. CHAUMET,<fnref n="1" />,* DANIEL SENTENAC,<fnref n="2" /> GUILLAUME MAIRE,<fnref n="1" /> MD RASEDUJJAMAN,<fnref n="1" /> TING ZHANG,<fnref n="3" /> ID AND ANNE SENTENAC<fnref n="1" />

Received 4 June 2021; revised 6 October 2021; accepted 25 October 2021; posted 28 October 2021; published 22 November 2021

In this tutorial, we introduce a solver of monochromatic Maxwell equations made freely available at [https://www.fresnel.fr/perso/chaumet/ifdda.html](https://www.fresnel.fr/perso/chaumet/ifdda.html), based on the volume moment method. The Institut Fresnel Discrete Dipole Approximation or Idiot-Friendly Discrete Dipole Approximation (IFDDA) calculates the diffracted field, the optical forces, and the image through a microscope of any three-dimensional inhomogeneous object, possibly anisotropic, placed in a stratified medium. In this method, only the object is meshed so the required memory space is kept to a minimum. We describe the principle and the potentialities of IFDDA and present comparisons with Mie theory and experimental data to assess the accuracy of the method. In addition, we provide a user guide for first steps with the solver. We hope that you will use and enjoy this numerical tool! © 2021 Optical Society of America

https://doi.org/10.1364/JOSAA.432685

## 1. INTRODUCTION

There are numerous methods that enable the rigorous study of the diffraction of an electromagnetic wave by an object of arbitrary form and relative permittivity. One quotes, for example the finite-difference time-domain method (FDTD), the finite-element method (FEM), the multiple multipole method (MMP), the surface integral equation method, or the volume integral method, also named the method of moment (MoM). The reader is referred to the article by F. M. Kahnert for the advantages and drawbacks of the most common methods [1].

In this paper, we focus on the method of moments in which the monochromatic Maxwell equations are cast into a volume integral equation and, more precisely, on its discretized version known as the discrete dipole approximation (DDA). This method was proposed by Purcell and Pennypacker [2] in 1973 to study the scattering and absorption of light by nonspherical dielectric grains. It has several advantages, especially when one is interested in free-space scattering of objects placed in a stratified medium. First, it is applicable to arbitrarily shaped, inhomogeneous, anisotropic, or metallic objects [3]. Second, the outgoing wave condition and the boundary conditions at the interfaces of the stratified medium are automatically satisfied so that the computation of the electric field is restricted to the

volume of the scatterer only. The number of unknowns and the memory requirement are thus usually smaller than that of FDTD or FEM, but the system of linear equations to be solved is dense. The interest and potentialities of DDA are well described in Refs. [4,5], and comparison with other numerical methods can be found in Refs. [6,7].

Several DDA codes are already available to the public. The first DDA solver, proposed by Draine and Flatau, simulates the light scattering and absorption by isolated or arrays of irregular particles placed in a homogeneous background [8,9]. A second code, developed by Yurkin and Hoekstra, simulates the light scattering by an object near a plane substrate [10]. DDA in homogeneous space is also available as a toolbox of MATLAB [11].

Here, we propose two codes that extend further the potentialities of DDA. IFDDAM (for Institut Fresnel DDA adapted to Multilayers) is able to simulate the interaction between any arbitrary incident field and an inhomogeneous, possibly anisotropic, object placed in a multilayer. IFDDA, the Institut Fresnel DDA or Idiot-Friendly DDA, for its part, is specifically optimized to deal with objects placed in a homogeneous medium. Both codes provide scattering and absorption cross sections, far-field and near-field maps, and the images obtained by several

<footnote n="1">¹Institut Fresnel, Aix Marseille Univ, CNRS, Centrale Marseille, Marseille, France</footnote>

<footnote n="2">²European Gravitational Observatory, 56021 Cascina (PI), Italy</footnote>

<footnote n="3">³College of Information Science and Electronic Engineering (ISEE), Zhejiang University, Zhejiang Provincial Key Laboratory of Information Processing, Communication and Networking (IPCAN), Hangzhou 310027, China</footnote>

<footnote>*Corresponding author: patrick.chaumet@fresnel.fr</footnote>
optical microscopes (bright field, dark field, phase, holography). In addition, IFDDA provides the optical forces and torques acting on the object. IFDDA is coded in FORTRAN with a C++ ergonomic interface and Qt displays. Its drop-down menus and preimplemented configurations enable the user to launch a calculation easily and get a physical insight into her/his issue rapidly.

## 2. DESCRIPTION OF THE MAXWELL EQUATIONS SOLVER

In this section, we describe the main steps of the volume integral method, or equivalently of the discretized dipole approximation that is used to solve Maxwell equations in IFDDAM.

### A. Configuration

Hereafter, we consider the diffraction of monochromatic waves of wavenumber $k_0 = 2\pi/\lambda$, where $\lambda$ is the wavelength by an object located in a homogeneous background; see Fig. 1(a) or inside a stratified medium, Fig. 1(b). We depict the reference stratified medium by its possibly complex scalar relative permittivity $\epsilon_{\text{ref}}(z)$, where $(x, y, z)$ is an orthonormal basis. Note that our time convention for the monochromatic waves is $\exp(-i\omega t)$ so that the permittivity imaginary part should always be positive.

The substrate corresponds to the semi-infinite medium when $z \to -\infty$. Its wavenumber is noted $k^- = n^-k_0$, where $n^- = \sqrt{\epsilon_{\text{ref}}(z \to -\infty)}$. The superstrate corresponds to the semi-infinite medium when $z \to +\infty$. Its wavenumber is noted $k^+ = n^+k_0$ with $n^+ = \sqrt{\epsilon_{\text{ref}}(z \to +\infty)}$. An inhomogeneous object of support $\Omega$ is introduced inside this reference medium, so that the relative permittivity of the whole system reads $\vec{\epsilon}(\mathbf{r}) = \epsilon_{\text{ref}}(z) + 4\pi\vec{\chi}(\mathbf{r})$, where $\vec{\chi}(\mathbf{r})$ is equal to zero outside $\Omega$. We also assume that both the stratified medium and the object under study are nonmagnetic, i.e., $\mu = \mu_0$ everywhere. Note

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** (a) Object in a homogeneous background: $\epsilon_{\text{ref}} = 1$; (b) object embedded in a multilayer; the reference permittivity $\epsilon_{\text{ref}}$ depends only on $z$. The multilayer comprises $L$ interfaces indexed from 1, corresponding to the lowest $z$, to $L$ corresponding to the highest $z$, and $L+1$ different media, indexed from 0, the substrate to $L$, the superstrate. In the preimplemented configurations, when the illumination comes from a source placed in far field (plane waves, Gaussian beams), the latter is always placed at $z = -\infty$, i.e., the beams always propagate toward positive $z$.</figcaption></figure>

that DDA has been extended to magnetic objects, [12,13] and to nonlinear susceptibility [14], but these options are outside the scope of this tutorial. To ease the exposition of the method, we consider that the electromagnetic wave illuminating the system stems from a source current $\mathbf{S}$. This source may generate plane waves, Gaussian beams, [15] or antenna-like fields in preimplemented configurations. When the source creating the incident field is in far field (to form plane wave illuminations), it is always placed in the substrate of the multilayer. In this case, $n^-$ should be real positive.

### B. Description of the Volume Integral Method

From Maxwell equations, the total electric field at $\mathbf{r} \in \mathbb{R}^3$ satisfies [16]

$$ \nabla \times \nabla \times \mathbf{E}(\mathbf{r}) - \epsilon_{\text{ref}}(z) k_0^2 \mathbf{E}(\mathbf{r}) = \mathbf{S} + 4\pi k_0^2 \vec{\chi}(\mathbf{r}) \mathbf{E}(\mathbf{r}). \quad (1) $$

The total field $\mathbf{E}(\mathbf{r})$ can be written as the sum of the reference field, $\mathbf{E}_{\text{ref}}(\mathbf{r})$, i.e., the field that would exist without the object, which verifies the homogeneous equation

$$ \nabla \times \nabla \times \mathbf{E}_{\text{ref}}(\mathbf{r}) - \epsilon_{\text{ref}}(z) k_0^2 \mathbf{E}_{\text{ref}}(\mathbf{r}) = \mathbf{S}, \quad (2) $$

and a diffracted field $\mathbf{E}_d(\mathbf{r}) = \mathbf{E}(\mathbf{r}) - \mathbf{E}_{\text{ref}}(\mathbf{r})$, which satisfies the outgoing wave boundary condition. To calculate the total field, we introduce the Green tensor of the reference medium $\vec{\mathbf{G}}$, solution of

$$ \nabla \times \nabla \times \vec{\mathbf{G}}(\mathbf{r}, \mathbf{r}') - \epsilon_{\text{ref}}(z) k_0^2 \vec{\mathbf{G}}(\mathbf{r}) = 4\pi k_0^2 \mathbf{I} \delta(\mathbf{r} - \mathbf{r}'), \quad (3) $$

where $\mathbf{I}$ denotes the unit tensor that satisfies outgoing boundary conditions. The expression of the Green tensor of a stratified medium can be found in Ref. [17]. The total field is then the solution of the self-consistent volume integral equation,

$$ \mathbf{E}(\mathbf{r}) = \mathbf{E}_{\text{ref}}(\mathbf{r}) + \int_{\Omega} \vec{\mathbf{G}}(\mathbf{r}, \mathbf{r}') \vec{\chi}(\mathbf{r}') \mathbf{E}(\mathbf{r}') \mathrm{d}\mathbf{r}'. \quad (4) $$

Once the field in $\Omega$ is known, it can be calculated everywhere.

### C. Estimating the Field inside the Object

To solve Eq. (4), we approximate the electric field and permittivity inside $\Omega$ by step-wise functions that are constant over small cubic subunits of side $d$ [18–22]. The choice of $d$ is crucial for the accuracy of the results, and it should be adapted to the spatial behavior of the field in the sample (exponential decay in conductive or absorptive materials, oscillation in dielectric material). As a rule of thumb, $d$ should be smaller than $\lambda/(2\pi|n|)$ [8,21,23]. Once discretized, the volume integral equation, restricted to $\mathbf{r} \in \Omega$, can be transformed into the linear system,

$$ \mathbf{E}(\mathbf{r}_i) = \mathbf{E}_{\text{ref}}(\mathbf{r}_i) + \sum_{j=1}^{N} \vec{\mathbf{g}}(\mathbf{r}_i, \mathbf{r}_j) \vec{\chi}(\mathbf{r}_j) \mathbf{E}(\mathbf{r}_j), \quad (5) $$

with $i = 1, \dots, N$, where $N$ denotes the number of nodes of the cubic mesh of $\Omega$, $\mathbf{E}(\mathbf{r}_i)$ the macroscopic field of the subunit $i$, and $\vec{\mathbf{g}}(\mathbf{r}_i, \mathbf{r}_j)$ the integral of $\vec{\mathbf{G}}(\mathbf{r}_i, \mathbf{r})$ over the subunit centered about $\mathbf{r}_j$, i.e., $\vec{\mathbf{g}}(\mathbf{r}_i, \mathbf{r}_j) = \int_{\Omega_j} \vec{\mathbf{G}}(\mathbf{r}_i, \mathbf{r}) \mathrm{d}\mathbf{r}$ where
$\Omega_j$ is the volume of the subunit $j$. For $\mathbf{r}_i \neq \mathbf{r}_j$ the integration of the Green tensor over the cubic subunit is usually approximated, assuming $\vec{\mathbf{G}}$ to be constant over the subunit, by $\vec{\mathbf{g}}(\mathbf{r}_i, \mathbf{r}_j) \approx d^3\vec{\mathbf{G}}(\mathbf{r}_i, \mathbf{r}_j)$. The integration of the diagonal term, $\mathbf{r}_i = \mathbf{r}_j$, is more delicate due to the singular behavior of the Green tensor at the origin; see Refs. [21,24,25]. One generally uses $\vec{\mathbf{g}}(\mathbf{r}_i, \mathbf{r}_i) = -\frac{4\pi}{3\varepsilon_{\text{ref}}} \vec{\mathbf{I}} + \text{PV}[\int_{\Omega_j} \vec{\mathbf{G}}(\mathbf{r}_i, \mathbf{r})d\mathbf{r}]$, where PV is the integration per principal value obtained by removing a small sphere at the origin calculated using different approximations following Ref. [25] and annexes 10 and 9.

The set of linear equations, Eq. (5), can be written symbolically as

$$ \mathbf{E} = \mathbf{E}_{\text{ref}} + \mathbf{A}\mathbf{D}_{\chi}\mathbf{E}, \qquad (6) $$

where $\mathbf{E}$ and $\mathbf{E}_{\text{ref}}$ are 3N vectors representing the unknown field inside $\Omega$ and the reference field, respectively. $\mathbf{A}$ is a $3N \times 3N$ matrix that contains the Green tensor, and $\mathbf{D}_{\chi}$ is a diagonal matrix of size $3N \times 3N$ that contains the tensor of permittivity contrast. The linear system is solved iteratively using a conjugate gradient method; see Ref. [22]. To fasten the calculation, all the matrix vector products appearing in the iterative method are performed using three-dimensional (3D) fast Fourier transforms (FFTs) for the homogeneous space [26] and using two-dimensional (2D) FFT [27] for the multilayer case (note that 3D FFT could also be implemented if the object is in the substrate or superstrate of a multilayer [28]). To this aim, the linear system is solved within an 'FFT box' $\Omega_{\text{FFT}}$, defined as the smallest box enclosing $\Omega$. In $\Omega_{\text{FFT}}$, the permittivity contrast of the 'useless' subunits are put to zero. Hereafter, the discretization step $d$ is defined from the number of subunits, $N_j$, which is set by the user, along the longest dimension of $\Omega_{\text{FFT}}$.

Once the field inside $\Omega$ is known, the near field in a domain surrounding the object is computed directly using the discretized version of Eq. (4); see Ref. [29]. On the other hand, the far field is obtained in a faster way by taking advantage of the far-field expression of the Green tensor, as seen in the following.

## D. Calculating the Field Far from the Object

In this paragraph, we estimate the diffracted field at the observation point $\mathbf{r}$ in the far field of the object (such that $r \gg r'$ and $r \gg r'^2/\lambda$ for any $r' \in \Omega$). The far-field Green tensor can be approximated by $\vec{\mathbf{G}}(\mathbf{r}, \mathbf{r}') \approx \frac{e^{ik^\pm r}}{r} \vec{\mathbf{G}}_{\text{ff}}(\mathbf{k}^\pm, \mathbf{r}')$, where the wave vector $\mathbf{k}^\pm$ is defined by $\mathbf{k}^\pm = k^\pm \hat{\mathbf{r}}$ with $k^\pm = n^\pm k_0$ with $\hat{\mathbf{r}} = \mathbf{r}/r$ [16], and the superscript (+) indicates that $\mathbf{r}$ is in the superstrate, whereas the superscript (-) indicates that $\mathbf{r}$ is in the substrate. Note that if $n^+ $ has an imaginary part, the diffracted far field is null in the superstrate.

The far field diffracted by the object in the direction $\mathbf{k}^\pm$ can be estimated with the integral

$$ \mathbf{E}_{\text{d}}(\mathbf{r}) = \frac{e^{ik^\pm r}}{r} \int_{\Omega} \vec{\mathbf{G}}_{\text{ff}}(\mathbf{k}^\pm, \mathbf{r}') \vec{\chi}(\mathbf{r}') \mathbf{E}(\mathbf{r}') d\mathbf{r}'. \qquad (7) $$

This technique is appropriate if only a few directions of observation are considered. Indeed, the integral in Eq. (7) is time-consuming, especially if the number of cubic subunits

forming the object, $N$, is large [30]. If the far field is to be calculated along many directions of observations, it is preferable to use another approach using FFTs [30].

Due to the translational invariance of the reference geometry (vacuum or multilayer) in the transverse $(x, y)$ plane, the far-field Green tensor satisfies $\vec{\mathbf{G}}_{\text{ff}}(\mathbf{k}^\pm, \mathbf{r}') = \vec{\mathbf{M}}^\pm(\mathbf{k}_\parallel, \mathbf{z}')e^{-ik_\parallel \cdot \mathbf{r}'}$, where $\mathbf{k}_\parallel$ is the projection of $\mathbf{k}^\pm$ on the transverse plane and $\mathbf{r}_\parallel$ is the projection of $\mathbf{r}$ on the transverse plane. Then, the diffracted far field can be estimated as $\mathbf{E}_d(\mathbf{r}) = \frac{e^{ik^\pm r}}{r} \mathbf{e}_d^\pm(\mathbf{k}_\parallel)$, where

$$ \mathbf{e}_d^\pm(\mathbf{k}_\parallel) = \int \vec{\mathbf{M}}^\pm(\mathbf{k}_\parallel, \mathbf{z}') \mathcal{F}_{2D}[\vec{\chi}(\mathbf{r}_\parallel, \mathbf{z}') \mathbf{E}(\mathbf{r}_\parallel, \mathbf{z}')](\mathbf{k}_\parallel) d\mathbf{z}', \quad (8) $$

and $\mathcal{F}_{2D}[f](\mathbf{k}_\parallel) = \int f(\mathbf{r}_\parallel) \exp(-i\mathbf{k}_\parallel \cdot \mathbf{r}_\parallel) d\mathbf{r}_\parallel$ denotes the 2D Fourier transform. Numerically, the estimation of $\mathbf{e}_d^\pm(\mathbf{k}_\parallel)$ is performed using 2D FFTs, and a sum over the $N_z$ layers along $z$ forming the object,

$$ \mathbf{e}_d^\pm(\mathbf{k}_\parallel) \approx d^3 \sum_{k=1}^{N_z} \vec{\mathbf{M}}^\pm(\mathbf{k}_\parallel, \mathbf{z}_k') \text{FFT}_{2D}[\vec{\chi}(\mathbf{r}_\parallel, \mathbf{z}_k') \mathbf{E}(\mathbf{r}_\parallel, \mathbf{z}_k')]. \quad (9) $$

The calculation of Eq. (8) is much faster than that of the integral in Eq. (7). Its counterpart is that the far field is estimated only along the specific directions of observation, satisfying $\mathbf{k}_\parallel = (i\Delta k, j\Delta k)$ for $i$ and $j = -K/2, \dots, K/2 - 1$ with $\Delta k = 2\pi/(Kd)$, and $K$ the even number of points in the FFT.

Once the far field is estimated, the time averaged differential cross section can be expressed as

$$ \frac{dC_{\text{ext}}^{\pm}}{d\Omega} = \frac{1}{2} c \varepsilon_0 n^{\pm} |\mathbf{e}_{\text{d}}^{\pm}(\mathbf{k}_{\parallel})|^2. \qquad (10) $$

If the object is in vacuum and illuminated by a plane wave with amplitude $E_{\text{ref}}$, the extinction, scattering, and absorption cross sections and the asymmetrical parameter are evaluated as

$$ C_{\text{ext}} = \int \frac{dC_{\text{ext}}}{d\Omega} d\Omega = \frac{k_0}{E_{\text{ref}}^2} \sum_{j=1}^{N} \text{Im} [\mathbf{E}_{\text{ref}}^*(\mathbf{r}_j) \cdot \mathbf{P}(\mathbf{r}_j)], \quad (11) $$

$$ C_{\text{abs}} = \frac{k_0}{E_{\text{ref}}^2} \sum_{j=1}^{N} \text{Im} [\mathbf{E}^*(\mathbf{r}_j) \cdot \mathbf{P}(\mathbf{r}_j)], \quad (12) $$

$$ C_{\text{sca}} = \frac{k_0^4}{E_{\text{ref}}^2} \int \left\| \sum_{j=1}^{N} [\mathbf{P}(\mathbf{r}_j) - \mathbf{n}(\mathbf{n} \cdot \mathbf{P}(\mathbf{r}_j))] e^{-ik_0 \cdot \mathbf{n} \cdot \mathbf{r}_j} \right\|^2 d\Omega, \quad (13) $$

$$ g = \frac{k_0^3}{C_{\text{sca}} E_{\text{ref}}^2} \times \int (\mathbf{n} \cdot \mathbf{k}_0) \left\| \sum_{j=1}^{N} [\mathbf{P}(\mathbf{r}_j) - \mathbf{n}(\mathbf{n} \cdot \mathbf{P}(\mathbf{r}_j))] e^{-ik_0 \cdot \mathbf{n} \cdot \mathbf{r}_j} \right\|^2 d\Omega, \quad (14) $$

with $\mathbf{P}(\mathbf{r}_j) = d^3\vec{\chi}(\mathbf{r}_j)\mathbf{E}(\mathbf{r}_j)$. The scattering cross section can be estimated with $C_{\text{sca}} = C_{\text{ext}} - C_{\text{abs}}$ or with the integral Eq. (13)
[18], where the integration over dΩ is either done in spherical coordinates [18] or in Cartesian coordinates, using FFT. This last approach is faster, but its accuracy depends on the precision with which the disk $k_\parallel < k_0$ is described with the square meshing of the FFT. This issue is particularly pressing at grazing incident angle, when most of the energy is scattered for $k_\parallel$ close to $k_0$. Note that the number K of points of the FFT can be set by the user (in the advanced interface).

## E. Simulating the Image of an Object Given by a Microscope

IFDDA comes with a microscope simulator. We consider a microscope consisting of an objective of numerical aperture (NA) satisfying the Abbe sine condition [31] and a tube lens in a 4f configuration leading to a magnification factor M. The z axis normal to the multilayer coincides with the optical axis. By convention, the object focal point is placed at $\mathbf{r} = \mathbf{r}_o = (0, 0, z_o)$, while the image focal plane is located at $z = z_i$. Recalling that the illuminating beam always propagates toward the positive z, the transmission configuration corresponds to $z_i > z_o$, Fig. 2(a), while the reflection configuration corresponds to $z_i < z_o$, Fig. 2(b). We assume that the medium between the tube lens and camera is always air. On the other hand, since the sample can be placed in a multilayer, the microscope simulator can account for any index mismatch between the objective oil, the coverslip, and the immersion medium.

The image given by the microscope is computed following the procedure detailed in Ref. [32]. Here, we briefly describe the main computational steps. The field in the substrate (–) or in the superstrate (+) propagating towards the objective in the reflection or transmission configuration, respectively, can be written as a sum of plane waves

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Sketch of the microscope. Left column, in transmission configuration, right column, in reflection configuration. $z_o$ denotes the position of the object focal plane and $z_i$ the position of the image focal plane. $\mathbf{k}$ is the wave vector of a plane wave diffracted by the object and $\mathbf{k}'$ is the wave vector of the same wave after transmission by the lenses. Note that the objective immersion medium (oil, water, or air) corresponds to the superstrate in the transmission configuration and to the substrate in reflection configuration.</figcaption></figure>

$$ \mathbf{E}_{\text{obj}}(\mathbf{r}) = \int \mathbf{e}_{\text{pupil}}^{\pm}(\mathbf{k}_{\parallel}) e^{i\mathbf{k}^{\pm} \cdot (\mathbf{r} - \mathbf{r}_o)} d\mathbf{k}_{\parallel}, \quad (15) $$

where

$$ \mathbf{e}_{\text{pupil}}^{\pm}(\mathbf{k}_{\parallel}) = \frac{\mathbf{e}_{\text{d}}^{\pm}(\mathbf{k}_{\parallel})}{-2i\pi|k_z|} e^{i\mathbf{k}_{\parallel}^{\pm}z_o} + \mathbf{e}_{\text{ref}}^{\pm}(\mathbf{k}_{\parallel}), \quad (16) $$

with $\mathbf{k}^{\pm} = \mathbf{k}_{\parallel} + k_z^{\pm} \hat{\mathbf{z}}$ and $k_z^{\pm} = \pm\sqrt{(n^2 k_0)^2 - k_{\parallel}^2}$ depending on the microscope configuration (transmission with the sign + and reflection with the sign –), $\mathbf{e}_{\text{d}}^{\pm}$ is given by Eq. (8), and $\mathbf{e}_{\text{ref}}^{\pm}$ represents the plane wave decomposition of the reference field following Eq. (15) with (–) indicating the reflected beam (propagating towards negative z) in the substrate and (+) indicating the transmitted beam in the superstrate. At the pupil or Fourier plane of the microscope, one has access to $\mathbf{e}_{\text{pupil}}^{\pm}(\mathbf{k}_{\parallel})$.

If $k_\parallel < k_0$ NA, the microscope transforms the plane wave with wave vector $\mathbf{k}$ in the object space into a plane wave with wave vector $\mathbf{k}'$ in the image space, with

$$ \mathbf{k}' = (\mathbf{k}'_{\parallel}, \pm k'_z) = \left( -\frac{\mathbf{k}_{\parallel}}{M}, \pm\sqrt{k_0^2 - \frac{\mathbf{k}_{\parallel}^2}{M^2}} \right). \quad (17) $$

The transformation of $\mathbf{k}$ into $\mathbf{k}'$ preserves the component of the electric field normal to the ($\mathbf{k}, \hat{\mathbf{z}}$) plane, (TE), and rotates the TM component (electric field in the ($\mathbf{k}, \hat{\mathbf{z}}$) plane). At the image focal plane, the field reads

$$ \mathbf{E}_{\text{ima}}(\mathbf{k}_{\parallel}, z_i) = \frac{1}{M} \int \sqrt{\frac{k_z}{k'_z}} \tilde{h}(\mathbf{k}_{\parallel}) \mathbf{R}(\mathbf{k}_{\parallel}) \mathbf{e}_{\text{pupil}}^{\pm}(\mathbf{k}_{\parallel}) e^{i\mathbf{k}'_{\parallel} \cdot \mathbf{r}_{\parallel}} d\mathbf{k}'_{\parallel}, \quad (18) $$

where $\tilde{h}(\mathbf{k}_{\parallel})$ is the pupil function, which acts as a low-pass filter, $\tilde{h}(\mathbf{k}_{\parallel}) = 1$ for $k_\parallel < k_0$ NA and 0 elsewhere. $\mathbf{R}(\mathbf{k}_{\parallel})$ is the rotation matrix defined as

$$ \mathbf{R}(\mathbf{k}_{\parallel}) = \begin{pmatrix} u_x^2(1 - \cos\theta) + \cos\theta & u_x u_y(1 - \cos\theta) & u_y \sin\theta \\ u_x u_y(1 - \cos\theta) & u_y^2(1 - \cos\theta) + \cos\theta & -u_x \sin\theta \\ -u_y \sin\theta & u_x \sin\theta & \cos\theta \end{pmatrix}, \quad (19) $$

where $\mathbf{u} = \frac{\hat{\mathbf{k}} \times \hat{\mathbf{z}}}{|\hat{\mathbf{k}} \times \hat{\mathbf{z}}|}$ is the rotation axis with $\hat{\mathbf{k}} = \mathbf{k}/k$. Note that $\mathbf{u}$ has no component along the z direction. $\theta$ is defined as $\cos\theta = \hat{\mathbf{k}} \cdot \hat{\mathbf{k}}'$ and $\sin\theta = \|\hat{\mathbf{k}} \times \hat{\mathbf{k}}'\|$. The factor $\frac{1}{M}\sqrt{\frac{k_z^2}{k'_z}}$ ensures the energy conservation; see appendix of Ref. [32].

## F. Evaluating the Transmissivity, Reflectivity, and Absorptivity

IFDDA allows one to estimate the reflectivity $\rho$, transmissivity $\tau$, and absorptivity $1 - \rho - \tau$ of the perturbed multilayered system [33]. These quantities are calculated only when the incident field stems from far-field sources (located in the substrate) so that it can be cast as a sum of plane waves propagating toward positive z. At $\mathbf{r}$ in the substrate, the incident field reads

$$ \mathbf{E}_{\text{inc}}(\mathbf{r}) = \int \mathbf{e}_{\text{inc}}(\mathbf{k}_{\parallel}) e^{i\mathbf{k} \cdot \mathbf{r}} d\mathbf{k}_{\parallel}, \quad (20) $$
with $k = k_{\parallel} + k_z \hat{z}$, $k_z = \sqrt{(n^{-}k_0)^2 - k_{\parallel}^2}$. The reflectivity and the transmissivity, $\rho$ and $\tau$, respectively, are estimated from the diffracted field, Eq. (8), as

$$ \tau = \frac{\int_{k_0 n^+} |\mathbf{e}_{\text{pupil}}^{+}(\mathbf{k}_{\parallel})|^2 d\mathbf{k}_{\parallel}}{\int_{k_0 n^-} |\mathbf{e}_{\text{inc}}(\mathbf{k}_{\parallel})|^2 d\mathbf{k}_{\parallel}} \quad \text{for } k_z > 0, \quad (21) $$

$$ \rho = \frac{\int_{k_0 n^-} |\mathbf{e}_{\text{pupil}}^{-}(\mathbf{k}_{\parallel})|^2 d\mathbf{k}_{\parallel}}{\int_{k_0 n^+} |\mathbf{e}_{\text{inc}}(\mathbf{k}_{\parallel})|^2 d\mathbf{k}_{\parallel}} \quad \text{for } k_z < 0. \quad (22) $$

When the sample and multilayer are non absorbing, one should check that $\rho + \tau$ is close to 1. Note that when the multilayer supports guided waves, the energy conservation cannot be checked, as the energy propagating in the guided wave is not computed; see Ref. [7].

### 3. COMPARISON OF IFDDA WITH MIE THEORY

The ability of IFDDAM to deal with complex configurations was checked with Maxwell solvers based on FDTD, finite elements, and the Fourier modal method. We considered a configuration where a grating coupler surrounded by Bragg mirrors, forming a $200\lambda \times 200\lambda \times \lambda/2$ structure, is deposited on a multilayer supporting guided modes. The agreement between IFDDAM and the other solvers at and outside the resonance conditions was remarkable; see Ref. [7].

In this section, we present comparisons of IFDDA with the analytical Mie theory [34], which is included in the code. This comparison is most useful to check if the discretization

<figure><img src="image_6.png" /><figcaption>Fig. 3. Relative errors between IFDDA and Mie simulations for a nonabsorbing sphere of radius $a = 10$ µm and three different discretization steps ($N$ is the number of subunits in $\Omega_{FFT}$) as a function of the sphere permittivity. The sphere is illuminated by a plane wave at $\lambda = 500$ nm. (a) Extinction cross section; (b) optical force; (c) asymmetry factor; and (d) the extinction cross section estimated with Mie theory and IFDDA with $N = 8 \times 10^6$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Same as Fig. 3, but the sphere permittivity is complex with $Re(\epsilon) = 2$ and the imaginary part is increased. (a) Extinction cross section; (b) absorbing cross section; (c) optical force; and (d) asymmetry factor.</figcaption></figure>

parameter $d$ has been well chosen. All the calculations of this section were performed using $\frac{2}{3}n_{\text{ref}}i k_0^3$ as the principal value of the Green function, i.e., the polarizability $\alpha_{RR}$; see Appendix D.

In Fig. 3, we considered a nonabsorbing sphere of radius $a = 10$ µm with increasing relative permittivity illuminated by a plane wave with $\lambda = 500$ nm. We compare the extinction cross section [34], the optical force and asymmetry factor to those provided by Mie theory for three different mesh sizes: $d = 100$ nm, $d = 66$ nm, and $d = 50$ nm, corresponding to $N = 8, 27$, and 64 millions of subunits for $\Omega_{FFT}$, respectively. As expected, the relative errors between IFDDA and Mie decrease with $d$ and with the sphere relative permittivity; see Fig. 3. The errors above 10% that are observed for large $d$ can be explained by the extreme sensitivity of the Mie resonances to the object

<figure><img src="image_8.png" /><figcaption>Fig. 5. (a) Relative errors between IFDDA and the Mie theory for the optical force and extinction cross section for an absorbing sphere of radius $a = \lambda/2$, $Im(\epsilon) = 1$ illuminated by a plane wave as a function of $Re(\epsilon)$ and two discretization steps, $N_l = 100$ and 200; (b) same as (a), but the sphere permittivity is set to 2 and the relative errors for the optical force, cross section, and asymmetrical parameter are plotted versus the discretization parameter $N_l$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Experimental setup of the holographic microscope in reflection configuration. HW, motorized half-wave plate; M, rotative mirror; BE, beam expander; BS, beam splitter; $L_1$, tube lens; $L_2$ and $L_3$, relay lenses; OL, objective lens.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Comparison of the experimental (left) and simulated (right) intensity recorded at the image plane of a microscope in reflection configuration; NA = 1.45 and magnification factor, $M$ = 200. The sample is a sphere of radius 2500 nm and relative permittivity of 2.5857, placed in water and deposited on a coverslip. The sphere is illuminated from the coverslip by a collimated beam making an angle $\theta$ = 30° with respect to the z axis. (a), (c) TE polarization; (b), (d) TM polarization.</figcaption></figure>

shape and the difficulty of reproducing the sphere curved surface with a cubic meshing. The oscillation of the cross-section error comes from a slight shift between the Mie and IFDDA cross-section curves; see Fig. 3(d).

In Fig. 4, the same test is conducted with an absorbing sphere, the real part of the relative permittivity being set to 2 and the imaginary part being increased. In this case, we observe that the accuracy of IFDDA is always better than 3%, the absorption damping the resonances.

To investigate further the limits of IFDDA, we plot in Fig. 5(a), the errors between IFDDA and Mie for a high refractive index sphere of radius $a = \lambda/2$. The real part of the sphere relative permittivity is increased from 1 to 30, with an imaginary part fixed to 1. We observe that the results obtained with

$N_l = 100$ and $N_l = 200$ are reasonably accurate (better estimations could be obtained by accounting for local-field corrections [35,36], but this option is only implemented in the code for spherical particles). Thus, DDA can handle samples with high permittivity, but it requires such a fine discretization that only small samples are computationally tractable. To give an insight on the convergence of IFDDA, the errors on $C_{ext}$, the total force and $g$ are plotted as a function of $N_l$ in Fig. 5(b) for a sphere of radius $a = 2\lambda$ and relative permittivity $\epsilon = 2$. A discussion on the accuracy of IFDDA with respect to the size of the object, the relative permittivity, and the mesh size can be found in Ref. [8].

## 4. COMPARISON OF IFDDAM WITH EXPERIMENTS

In this section, we compare the microscope images given by IFDDAM to those obtained experimentally with the microscope in reflection configuration described in Refs. [37,38] and illustrated in Fig. 6. The samples are placed on a coverslip of refractive index 1.5 that matches that of the oil objective of

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Same as Fig. 7, but the sample is now a resin star in air deposited on a coverslip and illuminated by a TE-polarized collimated beam (the incident field is directed along y) with $\theta = 68^\circ$ that is totally reflected at the glass-air interface and $z_o = 350$ nm. (a) Electron microscope image of the sample. The left (right) column corresponds to experimental (numerical) data. (b), (e) field intensity at the image plane; (c), (f) modulus of the y component of the diffracted field at the Fourier plane; (d), (g) phase of the y component of the diffracted field at the Fourier plane.</figcaption></figure>
NA = 1.45. They are illuminated by a collimated beam of wave-
length 475 nm, assimilated to a plane wave that comes from the
glass substrate and makes an angle θ with respect to the optical
axis. The plane of incidence is (x, z). An off-axis holographic
system allows us to measure the complex field at the image plane.

The first sample on the coverslip is a sphere of index 1.608 and diameter 5 µm surrounded by water (index 1.33). The microscope is tuned so that the object focal plane cuts the middle of the sphere at 2500 nm above the glass substrate. The magnification of the microscope, which should not be overlooked when dealing with polarized waves [32], is 200. Figure 7 displays the experimental and simulated intensities when the incident beam illuminates the sample at an angle θ = 30° with TE and TM polarizations. Experimentally, the position of the focal plane was estimated to z₀ = 2500 nm, with an accuracy of ±300 nm. We observe a good agreement for TE polarization with z₀ = 2500 nm and for TM polarization with z₀ = 2800 nm. The remaining differences may be explained by the uncertainty on the angle of incidence ±2°. The nonabsorbing glass sphere being a resonant object (as its optical volume is larger than one wavelength cube), a slight change in the illumination conditions may significantly affect the image pattern. To see the sensitivity of the image to the different parameters, see Appendix G.

The second sample on the coverslip consists of a resin star
made of 12 branches of width of 90 nm, length of 400 nm,
and height 160 nm of refractive index 1.5 in air. The inner
diameter of the star is 800 nm. A TE-polarized collimated beam
illuminates the sample at an angle of θ = 68° so that it is totally
reflected at the glass–air interface. The microscope is tuned
so that the object focal plane is placed in air at z₀ = 350 nm
above the substrate, and its magnification is 290. Figure 8 shows

the experimental and numerical images obtained at the image
focal plane and at the Fourier plane. We observe a very good
agreement between the simulation and the experiment.

**5. USING IFDDA AND IFDDAM IN PRACTICE**

The source code can be downloaded at Code 1, Ref. [39].

The package comes with a detailed user guide, which explains
how to use the interface, the configuration parameters, and
the data files. Two versions of the codes are available. IFDDA
is optimized for objects in free space and IFDDAM for objects
embedded in a multilayer. The free-space version is much
quicker than the multilayer counterpart and should be preferred
when appropriate. Both versions have a user-friendly interface
with a drop-down menu that has been designed so that nonspe-
cialists can easily perform a simulation and obtain meaningful
images (microscope images, near-field images, scattering cross
sections, among others), within a few minutes; see Fig. 9. All the
data can be displayed directly on the interface [40]. For exam-
ple, images of all the field components for all the possible cuts,
i.e., x, y, or z within the FFT box are available. The data are also
saved in ASCII or HDF5 format. A MATLAB code (ifdda.m) is
provided to read the HDF5 or ASCII data.

IFDDA comes with many preimplemented examples (see
user guide). The objects can be spheres, cuboids, ellipsoids,
and even random media. The incident field can be plane
waves, Gaussian beams [15], or antenna beams. For experts,
it is possible to provide custom incident fields and objects
using home-made ASCII files. For example, one can study the
scattering from inhomogeneous anisotropic objects embed-
ded in a waveguide and illuminated by a (provided) guided
mode or study the diffraction by a (provided) hundreds of

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Code interface. The drop-down menu on the left allows one to choose the object, illumination, and the wanted simulation (microscopy, optical forces, far-field scattering, etc.). The right side of the interface displays the data with color plots (near-field image, field at the image and Fourier planes of the microscope, force field, etc.).</figcaption></figure>
wavelength long, highly resonant, 3D subwavelength patterned structure [41].

Note that many numerical parameters and code options are set by default, but they can be modified in the advanced interface option panel.

A parallelized version using openmp (version >= 4.5) and the fastest Fourier transform in the west [42] is available for Linux systems. Another version (without the parallel calculation option) [43] is available for Windows.

# 6. CONCLUSION

Our ambition is that undergraduate students or researchers who are not specialists in electromagnetism try IFDDA or IFDDAM and its preimplemented examples (objects such as spheres or cuboids, incident fields such as Gaussian beams or plane waves) to get useful simulations (microscope images, optical forces) with a minimal effort.

Yet, we believe that these codes can also appeal to researchers in electromagnetism and microscopy, as they permit simulations of the light-matter interaction in very complex configurations.

# APPENDIX A: AUTHORS

P.C. developed most of the FORTRAN code and a part of the C++ code. He should be contacted first for any problem concerning the code. D.S. did the main part of the C++ code which manages the graphical interface and the makefile. G.M. and M.R. provided the experimental data. T.Z. optimized the computation of the diffracted field and A.S. helped on the far-field computation, the microscope simulator, and found many bugs.

# APPENDIX B: SOME REMARKS CONCERNING THE GREEN TENSORS

In the free-space version, IFDDA, the permittivity of the medium surrounding the object is 1; see Fig. 1(a). Of course, this choice does not limit the configurations to objects in vacuum. Indeed, the field of a monochromatic wave of wavelength $\lambda$ propagating through an object of relative permittivity contrast $\chi(\mathbf{r})$ in a lossless homogeneous background of permittivity $\epsilon_{\text{ref}}$ is the same as that of a monochromatic wave with wavelength $\lambda/\sqrt{\epsilon_{\text{ref}}} \text{ propagating through an object of relative permittivity } \vec{\epsilon}(\mathbf{r})/\epsilon_{\text{ref}} \text{ placed in vacuum. In this version, the Green tensor in homogeneous space is analytical; see Ref. [16]: its calculation is fast and easy.}$

In the multilayer version, IFDDAM, the key difficulty lies in the estimation of the Green tensor of the stratified reference medium. The calculation of the Green tensor in a stratified medium from the Weyl integral is made delicate by the possible presence of poles on the real axis corresponding to guided modes. To avoid them, the integration is performed in the complex plane following the approach proposed by [17]. By using an elliptical path starting at $k=0$ with the major semi-axis $k_{\text{max}} = k_0 n_{\text{max}} + k_0/2$ with $n_{\text{max}} = \max(\sqrt{\text{Re}(n)})$ and the minor semi-axis $k_{\text{min}} = 0.3 n_{\text{max}}$, the poles are avoided in all cases. However, the presence of exponentially growing functions in the T-matrix recursive algorithm yielding the Green tensor in the stratified medium limits the thickness of the multilayer to about $50\lambda$ (and even less if the layer is made of metal). This

issue could be avoided by replacing the T-matrix approach with a S-matrix technique [44].

Although efficient, the calculation of the Green tensor using this approach is time-consuming, especially if it has to be repeated for all the pairs $(\mathbf{r}, \mathbf{r}')$ with $(\mathbf{r}, \mathbf{r}') \in \Omega^2$. Due to the translational invariance of the geometry in the transverse $(x, y)$ plane, the Green tensor depending on $(|\mathbf{r}_\parallel - \mathbf{r}'_\parallel|, z, z')$, we interpolate $\mathbf{G}$ from a few rigorous calculations, $\mathbf{G}(qd/n_d, z, z')$ with $q = 1, \dots$, integer $(\sqrt{n_x^2 + n_y^2})$ and $n_d$ a natural number. Note that classical linear and polynomial interpolations cannot evaluate the Green tensor properly when $\|\mathbf{r}_\parallel - \mathbf{r}'_\parallel\| < \lambda$, as the fast decay of the evanescent waves is not accounted for accurately. We obtained better results using an interpolation with rational functions that are quotients of polynomials [41].

# APPENDIX C: MACROSCOPIC FIELD VERSUS LOCAL FIELD

In the near-field integral equation given by Eq. (5), the unknown field corresponds to the macroscopic field. Bearing in mind the singularity of the Green tensor, these equations can be rewritten with the local field as the unknown. More details about the relationship between the macroscopic and local-field approaches are given in Ref. [25]. If we neglect the principal value in the self-term of the Green tensor, i.e., $\mathbf{g}(\mathbf{r}_i, \mathbf{r}_j) \approx -\frac{4\pi}{3\epsilon_{\text{ref}}}$, the near-field equation becomes

$$ \mathbf{E}(\mathbf{r}_i) = \mathbf{E}_{\text{ref}}(\mathbf{r}_i) + \sum_{j=1, i \neq j}^{N} \vec{\mathbf{g}}(\mathbf{r}_i, \mathbf{r}_j) d^3 \vec{\chi}(\mathbf{r}_j) \mathbf{E}(\mathbf{r}_j) - \frac{4\pi \vec{\chi}(\mathbf{r}_i)}{3\epsilon_{\text{ref}}} \mathbf{E}(\mathbf{r}_i). \quad (C1) $$

Then, by gathering all the terms $\mathbf{E}(\mathbf{r}_i)$ to the left of the equation, the near-field equation can be written as

$$ \mathbf{E}_{\text{loc}}(\mathbf{r}_i) = \mathbf{E}_{\text{ref}}(\mathbf{r}_i) + \sum_{j=1, i \neq j}^{N} \vec{\mathbf{g}}(\mathbf{r}_i, \mathbf{r}_j) \vec{\alpha}_{\text{CM}}(\mathbf{r}_j) \mathbf{E}_{\text{loc}}(\mathbf{r}_j), \quad (C2) $$

where

$$ \mathbf{E}_{\text{loc}}(\mathbf{r}_i) = \left( \frac{\vec{\epsilon}(\mathbf{r}_j) + 2\epsilon_{\text{ref}}}{3\epsilon_{\text{ref}}} \right) \mathbf{E}(\mathbf{r}_i), \quad (C3) $$

denotes the local field, i.e., the field at the subunit position in the absence of the subunit itself and

$$ \vec{\alpha}_{\text{CM}}(\mathbf{r}_j) = \frac{3}{4\pi} d^3 \epsilon_{\text{ref}} (\vec{\epsilon}(\mathbf{r}_j) - \epsilon_{\text{ref}}) (\vec{\epsilon}(\mathbf{r}_j) + 2\epsilon_{\text{ref}})^{-1} \quad (C4) $$

denotes the polarizability of each subunit related to the Clausius-Mossotti (CM) relation [18,20]. It is sometimes convenient to use this expression, as $\mathbf{p} = \vec{\alpha} \mathbf{E}_{\text{loc}}$ represents the dipolar moment associated with the subunits. Note that the expression of the polarizability directly depends on the accuracy with which the self-term $\mathbf{g}(\mathbf{r}_i, \mathbf{r}_j)$ is estimated, with important consequences for the field calculation, as shown hereafter.

# APPENDIX D: POLARIZABILITY

The notion of polarizability of the subunit introduced in the previous paragraph has a long history. The first to have been used is known as the CM polarizability and stems from the
simplest approximation of the self-term $g(\mathbf{r}_i, \mathbf{r}_j) \approx -\frac{4\pi}{3\epsilon_{\text{ref}}}$ [2]; see Eq. (26). Unfortunately, this simple approximation does not ensure energy conservation. The latter is obtained without neglecting the principal value in the self-term, $g(\mathbf{r}_i, \mathbf{r}_j) = -\frac{4\pi}{3\epsilon_{\text{ref}}} + \text{PV}[\int_{\Omega_i} \vec{G}(\mathbf{r}_i, \mathbf{r})d\mathbf{r}]$. Assuming the cubic subunit to be a sphere of the same volume, we get $\text{PV}[\int_{\Omega_i} \vec{G}(\mathbf{r}_i, \mathbf{r})d\mathbf{r}] = \frac{2}{3}n_{\text{ref}}ik_0^3$, [21]. With this better approximation, the polarizability reads [18,25]

$$ \alpha_{RR} = \frac{\alpha_{CM}}{1 - \frac{2}{3}ik_0^3 n_{\text{ref}} \alpha_{CM}}. \qquad (\text{D1}) $$

Other expressions of the self-term and subsequently of the polarizability have been proposed in order to improve the precision of the DDA and take into account the nonpunctual character of the dipole. Among the best known, one quotes that of Goedecke and O'Brien [45]

$$ \alpha_{GB} = \frac{\alpha_{CM}}{1 - \frac{2}{3}ik_0^3 n_{\text{ref}}\alpha_{CM} - k_0^2\alpha_{CM}/a}, \qquad (\text{D2}) $$

that of Lakhtakia [25]

$$ \alpha_{LA} = \frac{\alpha_{CM}}{1 - 2\frac{\epsilon - \epsilon_{\text{ref}}}{\epsilon + 2\epsilon_{\text{ref}}}[(1 - ik_0 n_{\text{ref}} a)e^{ik_0 n_{\text{ref}} a} - 1]}, \qquad (\text{D3}) $$

with $\frac{4\pi}{3}a^3 = d^3$ and Draine and Goodman [20]

$$ \alpha_{LR} = \frac{\alpha_{CM}}{1 + \alpha_{CM} \left[ \frac{(b_1 + \epsilon b_2 / \epsilon_{\text{ref}} + \epsilon b_3 / \epsilon_{\text{ref}} S) k_0^2}{d} - \frac{2}{3} i n_{\text{ref}} k_0^3 \right]}, \qquad (\text{D4}) $$

with $b_1 = -1.891531$, $b_2 = 0.1618469$, $b_3 = -1.7700004$, and $S = 1/5$. These different expressions have been included in the code (in the advanced interface). Note that the cubic mesh can easily be replaced by a cuboid mesh by changing the expression of the polarizability [46].

## APPENDIX E: OPTICAL FORCES

The free-space version, IFDDA, is able to compute the optical forces and torques acting on the object. The force on each subunit is computed as [47]

$$ F_u(\mathbf{r}_i) = \frac{1}{2} \left( \sum_{v=1}^{3} p_v(\mathbf{r}_i) \frac{\partial E_v^*(\mathbf{r}_i)}{\partial u} \right), \qquad (\text{E1}) $$

where $u$ and $v$ stand for either $x$, $y$, or $z$ and $*$ denotes the complex conjugate. The net force is given by $\mathbf{F} = \sum_i \mathbf{F}(\mathbf{r}_i)$. To perform the computation of the field derivative efficiently, we use FFTs, as detailed in Ref. [48]. The density of optical torque is calculated as [49]

$$ \Gamma(\mathbf{r}_i) = \mathbf{r}_i^g \times \mathbf{F}(\mathbf{r}_i) + \frac{1}{2} \text{Re} \left[ \mathbf{p}(\mathbf{r}_i) \times \left[ \mathbf{p}(\mathbf{r}_i) / \alpha_{\text{CM}}(\mathbf{r}_i) \right]^* \right], \qquad (\text{E2}) $$

where $\mathbf{r}_i^g$ is the vector between $i$ and the center of mass of the object and the net torque is obtained through $\Gamma = \sum_i \Gamma(\mathbf{r}_i)$. The calculation of the optical forces in presence of a substrate [50] is possible, but this option has not yet been implemented in IFDDAM.

## APPENDIX F: APPROXIMATE METHODS

The solving of the set of linear equations [Eq. (5)] required to estimate the field inside $\Omega$ is time-consuming and, in some cases, some approximate solutions are sufficient [51]. IFDDA and IFDDAM come with different approximate methods for estimating the field inside the object.

When the object size is small compared to the wavelength and the permittivity contrast between the object and the surrounding medium is also small, one can assume that the field inside $\Omega$ is not perturbed by the object. In that case, known as Born or single-scattering approximation, the field inside $\Omega$ is given by the reference field, $\mathbf{E} \approx \mathbf{E}_{\text{ref}}$. Another implementation of the Born approximation (called renormalized Born approximation; see Ref. [52]), which is usually more precise, consists in approximating the local field by the reference field, $\mathbf{E}_{\text{loc}} \approx \mathbf{E}_{\text{ref}}$ so that $\mathbf{E} \approx \frac{3\epsilon_{\text{ref}}}{\epsilon + 2\epsilon_{\text{ref}}} \mathbf{E}_{\text{ref}}$. Double scattering can be introduced using the next term in the Born series, $\mathbf{E} \approx \mathbf{E}_{\text{Born},1} = \mathbf{E}_{\text{ref}} + \mathbf{A}\mathbf{D}_\chi \mathbf{E}_{\text{ref}}$, at the cost of one matrix vector product.

For weakly contrasted objects that are large compared to the wavelength, other approximate solutions are available

<figure><img src="image_13.png" /><figcaption>Fig. 10. Simulated intensities recorded at the image plane of a microscope for different positions $z_o$ of the focal plane with respect to the coverslip surface set at $z=0$. The microscope is in reflection configuration, NA = 1.45, and magnification factor, $M=200$. The sample is a sphere of radius 2500 nm and relative permittivity of 2.5857 deposited on the coverslip and surrounded by water. The sphere is illuminated from the coverslip by a TE- or TM-polarized collimated beam with $\theta = 30^\circ$ with respect to the $z$ axis. Left column, TE polarization; right column, TM polarization; first line $z_o = 2200$ nm; second line $z_o = 2500$ nm; third line $z_o = 2800$ nm.</figcaption></figure>
in the free-space version, IFDDA. The field inside $\Omega$ can be approximated using the Rytov model [53,54]. In this case, the component $\beta$ of the field inside the object is written as

$$E^{\beta}(\mathbf{r}_i) \approx E_{\text{ref}}^{\beta}(\mathbf{r}_i) e^{E_{\text{Born},1}^{\beta}(\mathbf{r}_i)/E_{\text{ref}}^{\beta}(\mathbf{r}_i)}, \quad (\text{F1})$$

where $\beta = x, y, z$. Note that Rytov does not account for any field depolarization and needs one matrix vector product.

For even larger objects with mainly forward scattering, the beam propagation method (BPM) can be an interesting solution. BPM estimates the electromagnetic field inside the object by alternatively evaluating the diffraction and refraction inside the object. It is important to note that BPM ignores reflections; for more details, see Refs. [55,56]. The BPM field in the object is calculated z-slice by z-slice following the recurrence relation

$$\mathbf{E}(\mathbf{r}_{\parallel}, z + d) = e^{ik_0 n(\mathbf{r}_{\parallel}, z+d)d} \times \mathcal{F}_{2D}^{-1}[\mathcal{F}_{2D}[\mathbf{E}(\mathbf{r}_{\parallel}, z)](\mathbf{k}_{\parallel})e^{-i(k_0 - k_z)d}], \quad (\text{F2})$$

where $\mathcal{F}_{2D}$ is the 2D Fourier transform in the $(x, y)$ plane. In all cases, once the field inside the object is known, the diffracted field is obtained using Eq. (8).

<figure><img src="image_14.png" /><figcaption>Fig. 11. Same as Fig. 10, but the focal plane is set at $z_o = 2500$ nm and the incident angle is changed. Left column, TE polarization; right column, TE polarization; first line, $\theta = 28^{\circ}$; second line $\theta = 30^{\circ}$; third line $\theta = 32^{\circ}$.</figcaption></figure>

## APPENDIX G: INFLUENCE OF THE FOCAL PLANE POSITION ON THE IMAGE

The image provided by the microscope in Fig. 10 depends on the position of the focal plane $z_o$ with respect to the coverslip surface (placed at $z=0$). Experimentally, $z_o$ is estimated to be $2500 \pm 300$ nm. In Fig. 10, we simulated the images for $z_o = 2200, 2500$, and $2800$ nm. (Note that the field inside the sphere being unchanged, these simulations could be done very quickly). The sensitivity of the image to $z_o$, especially in TM polarization, explains partly the difficulty to perfectly match the simulation with the experiment.

In Fig. 11, we plot the intensity at the image plane for different incident angles $\theta = 28^{\circ}, 30^{\circ},$ and $32^{\circ}$. The position of the focal plane is set at $z_o = 2500$ nm. We observe that the intensity pattern in TM polarization is highly sensitive to a small variation of the incident angle. This analysis points out the dependence of the image on various parameters and the interest of a rigorous microscope simulator for studying the image formation.

## APPENDIX H: DISPLAYING 3D DATA WITH IFDDA

IFDDA provides the field inside the FFT box ($\Omega_{\text{FFT}}$) enclosing the target (and even beyond this box if one wants to study the field on a larger domain). In Fig. 12, we display the modulus of the macroscopic field inside $\Omega_{\text{FFT}}$ for multiple $(x, y)$ cuts at different $z$. The sample and illumination are the same as in

<figure><img src="image_15.png" /><figcaption>Fig. 12. Modulus of the electromagnetic field on multiple $(x, y)$ cuts at different $z$ of the FFT box $\Omega_{\text{FFT}}$ surrounding the sample. The configuration is the same as that of Fig. 10. The sphere is illuminated by a TM-polarized beam of incidence $\theta = 30^{\circ}$. $z=0$ corresponds to the coverslip surface. The size of the images is $5 \times 5$ $\mu$m$^2$, and each image has its own color scale.</figcaption></figure>
Fig. 10. These images were made using MATLAB from the HDF5 files provided by IFDDA, but similar plots are readily available on the IFDDA interface.

**Disclosures.** The authors declare no conflicts of interest.

**Data Availability.** Data underlying the results presented in this paper have all been realized by the IFDDA and IFDDAM codes, which can be downloaded at Code 1, Ref. [39].

## REFERENCES

1. F. M. Kahnert, "Numerical methods in electromagnetic scattering theory," J. Quant. Spectrosc. Radiat. Transfer 79-80, 775-824 (2003).

2. E. M. Purcell and C. R. Pennypacker, "Scattering and absorption of light by nonspherical dielectric grains," Astrophys. J. 186, 705-714 (1973).

3. M. A. Yurkin, D. de Kanter, and A. G. Hoekstra, "Accuracy of the discrete dipole approximation for simulation of optical properties of gold nanoparticles," J. Nanophoton. 4, 1-15 (2010).

4. M. A. Yurkin and A. G. Hoekstra, "The discrete dipole approximation: an overview and recent developments," J. Quant. Spectrosc. Radiat. Transfer 106, 558-589 (2007).

5. M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, "The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength," J. Quant. Spectrosc. Radiat. Transfer 106, 546-557 (2007).

6. M. A. Yurkin, A. G. Hoekstra, R. S. Brock, and J. Q. Lu, "Systematic comparison of the discrete dipole approximation and the finite difference time domain method for large dielectric scatterers," Opt. Express 15, 17902-17911 (2007).

7. P. C. Chaumet, G. Demésy, O. Gauthier-Lafaye, A. Sentenac, E. Popov, and A.-L. Fehrembach, "Electromagnetic modeling of large subwavelength-patterned highly resonant structures," Opt. Lett. 41, 2358-2361 (2016).

8. B. T. Draine and P. J. Flatau, "Discrete-dipole approximation for scattering calculations," J. Opt. Soc. Am. A 11, 1491-1499 (1994).

9. B. T. Draine and P. J. Flatau, "Discrete-dipole approximation for periodic targets: theory and tests," J. Opt. Soc. Am. A 25, 2693-2703 (2008).

10. M. A. Yurkin and A. G. Hoekstra, "The discrete-dipole-approximation code ADDA: capabilities and known limitations," J. Quant. Spectrosc. Radiat. Transfer 112, 2234-2247 (2011).

11. V. L. Loke, M. P. Menguc, and T. A. Nieminen, "Discrete-dipole approximation with surface interaction: computational toolbox for MATLAB," J. Quant. Spectrosc. Radiat. Transfer 112, 1711-1725 (2011), Special Issue on Electromagnetic and Light Scattering by Nonspherical Particles XII.

12. P. C. Chaumet and A. Rahmani, "Coupled-dipole method for magnetic and negative refraction materials," J. Quant. Spectrosc. Radiat. Transfer 110, 22-29 (2009).

13. P. C. Chaumet and A. Rahmani, "Electromagnetic force and torque on magnetic and negative-index scatterers," Opt. Express 17, 2224-2234 (2009).

14. N. K. Balla, P. T. So, and C. J. Sheppard, "Coupled dipole model for nonlinear scattering," in Frontiers in Optics 2009/Laser Science XXVI/Fall 2009, OSA Optics & Photonics Technical Digest (Optical Society of America, 2009), paper LSWK5.

15. P. C. Chaumet, "Fully vectorial highly non paraxial beam close to the waist," J. Opt. Soc. Am. A 23, 3197-3202 (2006).

16. J. D. Jackson, Classical Electrodynamics, 2nd ed. (Wiley, 1975).

17. M. Paulus, P. Gay-Balmaz, and O. J. F. Martin, "Accurate and efficient computation of the Green's tensor for stratified media," Phys. Rev. E 62, 5797-5807 (2000).

18. B. T. Draine, "The discrete-dipole approximation and its application to interstellar graphite grains," Astrophys. J. 333, 848-872 (1988).

19. P. J. Flatau, G. L. Stephens, and B. T. Draine, "Light scattering by rectangular solids in the discrete-dipole approximation: a new algorithm exploiting the Block-Toeplitz structure," J. Opt. Soc. Am. A 7, 593-600 (1990).

20. B. T. Draine and J. Goodman, "Beyond Clausius-Mossotti: wave propagation on a polarizable point lattice and the discrete dipole approximation," Astrophys. J. 405, 685-697 (1993).

21. P. C. Chaumet, A. Sentenac, and A. Rahmani, "Coupled dipole method for scatterers with large permittivity," Phys. Rev. E 70, 036606 (2004).

22. P. C. Chaumet and A. Rahmani, "Efficient iterative solution of the discrete dipole approximation for magneto-dielectric scatterers," Opt. Lett. 34, 917-919 (2009).

23. E. Zubko, D. Petrov, Y. Grynkov, Y. Shkuratov, H. Okamoto, K. Muinonen, T. Nousiainen, H. Kimura, T. Yamamoto, and G. Videen, "Validity criteria of the discrete dipole approximation," Appl. Opt. 49, 1267-1279 (2010).

24. A. D. Yaghjian, "Electric dyadic Green's functions in the source region," Proc. IEEE 68, 248-263 (1980).

25. A. Lakhtakia, "Strong and weak forms of the method of moments and the coupled dipole method for scattering of time-harmonic electromagnetics fields," Int. J. Mod. Phys. C 3, 583-603 (1992).

26. J. J. Goodman and P. J. Flatau, "Application of fast-Fourier-transform techniques to the discrete-dipole approximation," Opt. Lett. 16, 1198-1200 (2002).

27. R. Schmehl, B. M. Nebeker, and E. D. Hirleman, "Discrete-dipole approximation for scattering by features on surfaces by means of a two-dimensional fast Fourier transform technique," J. Opt. Soc. Am. A 14, 3026-3036 (1997).

28. M. A. Yurkin and M. Huntemann, "Rigorous and fast discrete dipole approximation for particles near a plane interface," J. Phys. Chem. C 119, 29088-29094 (2015).

29. P. J. Flatau and B. T. Draine, "Fast near field calculations in the discrete dipole approximation for regular rectilinear grids," Opt. Express 20, 1247-1252 (2012).

30. P. C. Chaumet, T. Zhang, and A. Sentenac, "Fast far-field calculation in the discrete dipole approximation," J. Quant. Spectrosc. Radiat. Transfer 165, 88-92 (2015).

31. E. Abbe, "VII.—On the estimation of aperture in the microscope," J. R. Microsc. Soc. 1, 388-423 (1981).

32. S. Khadir, P. C. Chaumet, G. Baffou, and A. Sentenac, "Quantitative model of the image of a radiating dipole through a microscope," J. Opt. Soc. Am. A 36, 478-484 (2019).

33. W. Weinstein, "The reflectivity and transmissivity of multiple thin coatings," J. Opt. Soc. Am. B 37, 576-581 (1947).

34. J. A. Stratton, Electromagnetic Theory (McGraw-Hill, 1941).

35. A. Rahmani, P. C. Chaumet, and G. W. Bryant, "On the importance of local-field corrections for polarizable particles on a finite lattice: application to the discrete dipole approximation," Astrophys. J. 607, 873-878 (2004).

36. A. Sentenac, P. C. Chaumet, and G. Leuchs, "Total absorption of light by a nanoparticle: an electromagnetic sink in the optical regime," Opt. Lett. 38, 818-820 (2013).

37. G. Maire, Y. Ruan, T. Zhang, P. C. Chaumet, H. Giovannini, D. Sentenac, A. Talneau, K. Belkebir, and A. Sentenac, "High-resolution tomographic diffractive microscopy in reflection configuration," J. Opt. Soc. Am. A 30, 2133-2139 (2013).

38. G. Maire, H. Giovannini, A. Talneau, P. C. Chaumet, K. Belkebir, and A. Sentenac, "Phase imaging and synthetic aperture super-resolution via total internal reflection microscopy," Opt. Lett. 43, 2173-2176 (2018).

39. P. C. Chaumet, A. Sentenac, and D. Sentenac, "Idiot Friendly-Discrete Dipole Approximation IF-DDA(M): electromagnetic scattering in three dimension," Github (2019) [accessed November 16, 2021], https://www.fresnel.fr/perso/chaumet/ifdda.html.

40. J. Blanchette and M. Summerfield, C++ GUI Programming with Qt 4, 2nd ed. (Pearson Education, 2008).

41. P. C. Chaumet, A. Sentenac, and A.-L. Fehrembach, "Modelling of hundreds of wavelength long, highly resonant, 3D subwavelength patterned scattering structures," Opt. Quantum Electron. 49, 71 (2017).

42. M. Frigo and S. G. Johnson, "The design and implementation of FFTW3," Proc. IEEE 93, 216-231 (2005), Special Issue on Program Generation, Optimization, and Platform Adaptation.

43. R. C. Singleton, "On computing the fast Fourier transform," Commun. ACM 10, 647-654 (1967).
44. L. Li, "Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings," *J. Opt. Soc. Am. A* **13**, 1024-1035 (1996).

45. G. H. Goedecke and S. G. O'Brien, "Scattering by irregular inhomogeneous particles via the digitized Green's function algorithm," *Appl. Opt.* **27**, 2431-2438 (1988).

46. D. A. Smunev, P. C. Chaumet, and M. A. Yurkin, "Rectangular dipoles in the discrete dipole approximation," *J. Quant. Spectrosc. Radiat. Transfer* **156**, 67-79 (2015).

47. P. C. Chaumet and M. Nieto-Vesperinas, "Time-averaged total force on a dipolar sphere in an electromagnetic field," *Opt. Lett.* **25**, 1065-1067 (2000).

48. P. C. Chaumet, A. Rahmani, A. Sentenac, and G. W. Bryant, "Efficient computation of optical forces with the coupled dipole method," *Phys. Rev. E* **72**, 046708 (2005).

49. P. C. Chaumet and C. Billaudieu, "Coupled dipole method to compute optical torque: application to a micropropeller," *J. Appl. Phys.* **101**, 023106 (2007).

50. P. C. Chaumet, A. Rahmani, and M. Nieto-Vesperinas, "Photonic force spectroscopy on metallic and absorbing nanoparticles," *Phys. Rev. B* **71**, 045425 (2005).

51. P. C. Chaumet, A. Sentenac, and T. Zhang, "Reflection and transmission by large inhomogeneous media. Validity of Born, Rytov and beam propagation methods," *J. Quant. Spectrosc. Radiat. Transfer* **243**, 106816 (2020).

52. K. Belkebir, P. C. Chaumet, and A. Sentenac, "Superresolution in total internal reflection tomography," *J. Opt. Soc. Am. A* **22**, 1889-1897 (2005).

53. T. M. Habashy, R. W. Groom, and B. R. Spies, "Beyond the Born and Rytov approximations–a nonlinear approach to electromagnetic scattering," *J. Geophys. Res. [Solid Earth Planets]* **98**, 1759–1775 (1993).

54. R. Carminati, "Phase properties of the optical near field," *Phys. Rev. E* **55**, R4901-R4904 (1997).

55. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, "Optical tomographic image reconstruction based on beam propagation and sparse regularization," *IEEE Trans. Comput. Imaging* **2**, 59-70 (2016).

56. J. V. Roey, J. van der Donk, and P. E. Lagasse, "Beam-propagation method: analysis and assessment," *J. Opt. Soc. Am.* **71**, 803-810 (1981).