## Abstract

In this tutorial, we introduce a solver of monochromatic Maxwell equations made freely available at 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

## 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 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 ${\varepsilon _{{\rm 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 {{\varepsilon _{{\rm 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 {{\varepsilon _{{\rm 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 $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \varepsilon} (\textbf{r}) = {\varepsilon _{{\rm ref}}}(z) + 4\pi \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \chi} (\textbf{r})$, where $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \chi} (\textbf{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 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 $\textbf{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 $\textbf{r} \in {{\mathbb R}^3}$ satisfies [16]

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

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 $\textbf{r} \in \Omega$, can be transformed into the linear system,

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

where $\textbf{E}$ and ${\textbf{E}_{{\rm ref}}}$ are $3N$ vectors representing the unknown field inside $\Omega$ and the reference field, respectively. $\textbf{A}$ is a $3N \times 3N$ matrix that contains the Green tensor, and ${\textbf{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 _{{\rm FFT}}}$, defined as the smallest box enclosing $\Omega$. In ${\Omega _{{\rm 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_l}$, which is set by the user, along the longest dimension of ${\Omega _{{\rm 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 $\textbf{r}$ in the far field of the object (such that $r \gg r^\prime $ and $r \gg {r^{\prime 2}}/\lambda$ for any $r^\prime \in \Omega$). The far-field Green tensor can be approximated by $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over {\textbf{G}}}(\textbf{r},\textbf{r}^\prime) \approx \frac{{{e^{i{k^ \pm}r}}}}{r}{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over {\textbf{G}}}_{{\rm ff}}}({\textbf{k}^ \pm},\textbf{r}^\prime)$, where the wave vector ${\textbf{k}^ \pm}$ is defined by ${\textbf{k}^ \pm} = {k^ \pm}\hat{\textbf{r}}$ with ${k^ \pm} = {n^ \pm}{k_0}$ with $\hat{\textbf{r}} = \textbf{r}/r$ [16], and the superscript ($+$) indicates that $\textbf{r}$ is in the superstrate, whereas the superscript ($-$) indicates that $\textbf{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 ${\textbf{k}^ \pm}$ can be estimated with the integral

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 ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over {\textbf{G}}}_{{\rm ff}}}({\textbf{k}^ \pm},\textbf{r}^\prime) = {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over {\textbf{M}}}^ \pm}({\textbf{k}_\parallel},z^\prime){e^{- i{\textbf{k}_\parallel} \cdot {\textbf{r}_\parallel}^\prime}}$, where ${\textbf{k}_\parallel}$ is the projection of ${\textbf{k}^ \pm}$ on the transverse plane and ${\textbf{r}_\parallel}$ is the projection of $\textbf{r}$ on the transverse plane. Then, the diffracted far field can be estimated as ${\textbf{E}_{\rm d}}(\textbf{r}) = \frac{{{e^{i{k^ \pm}r}}}}{r}\textbf{e}_{\rm d}^ \pm ({\textbf{k}_\parallel})$, where

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 ${\textbf{k}_\parallel} = (i\Delta k,j\Delta k)$ for $i$ and $j = - K/2, \cdots, 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

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

#### 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 $\textbf{r} = {\textbf{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} \gt {z_o}$, Fig. 2(a), while the reflection configuration corresponds to ${z_i} \lt {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

If ${k_\parallel} \lt {k_0}{\rm NA}$, the microscope transforms the plane wave with wave vector $\textbf{k}$ in the object space into a plane wave with wave vector $\textbf{k}^\prime $ in the image space, with

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

#### 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 $\textbf{r}$ in the substrate, the incident field reads

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 parameter $d$ has been well chosen. All the calculations of this section were performed using $\frac{2}{3}{n_{{\rm ref}}}ik_0^3$ as the principal value of the Green function, i.e., the polarizability ${\alpha _{{\rm RR}}}$; see Appendix D.

In Fig. 3, we considered a nonabsorbing sphere of radius $a = 10\;\unicode{x00B5}{\rm m}$ with increasing relative permittivity illuminated by a plane wave with $\lambda = 500 \; {\rm 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 \; {\rm nm}$, $d = 66 \; {\rm nm}$, and $d = 50 \; {\rm nm}$, corresponding to $N = 8$, 27, and 64 millions of subunits for ${\Omega _{{\rm 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 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_{{\rm 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 $\varepsilon = 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 ${\rm NA} = 1.45$. They are illuminated by a collimated beam of wavelength 475 nm, assimilated to a plane wave that comes from the glass substrate and makes an angle $\theta$ 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 $\theta = {30^\circ}$ with TE and TM polarizations. Experimentally, the position of the focal plane was estimated to ${z_o} = 2500 \; {\rm nm}$, with an accuracy of $\pm300 \;{\rm nm}$. We observe a good agreement for TE polarization with ${z_o} = 2500 \; {\rm nm}$ and for TM polarization with ${z_o} = 2800 \; {\rm nm}$. The remaining differences may be explained by the uncertainty on the angle of incidence $\pm{2^ \circ}$. 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 $\theta = {68^\circ}$ 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_0} = 350 \; {\rm 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 nonspecialists 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 example, 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 embedded in a waveguide and illuminated by a (provided) guided mode or study the diffraction by a (provided) hundreds of 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 (${\rm version}\gt = {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 (\textbf{r})$ in a lossless homogeneous background of permittivity ${\varepsilon _{{\rm ref}}}$ is the same as that of a monochromatic wave with wavelength $\lambda /\sqrt {{\varepsilon _{{\rm ref}}}}$ propagating through an object of relative permittivity $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \varepsilon} (\textbf{r})/{\varepsilon _{{\rm ref}}}$ 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_{\max}} = {k_0}{n_{\max}} + {k_0}/2$ with ${n_{\max}} = {\max}({\sqrt {{\rm Re}(n)}})$ and the minor semi-axis ${k_{\min}} = 0.3{n_{\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 $(\textbf{r},\textbf{r}^\prime)$ with $(\textbf{r},\textbf{r}^\prime) \in {\Omega ^2}$. Due to the translational invariance of the geometry in the transverse $(x,y)$ plane, the Green tensor depending on $(|{\textbf{r}_\parallel} - {\textbf{r}^\prime _\parallel}|,z,z^\prime)$, we interpolate $\textbf{G}$ from a few rigorous calculations, $\textbf{G}(qd/{n_d},z,z^\prime)$ with $q = 1, \cdots ,{\rm 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 $\parallel {\textbf{r}_\parallel} - {\textbf{r}_\parallel^\prime} \parallel \lt \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., $\textbf{g}({\textbf{r}_i},{\textbf{r}_i}) \approx - \frac{{4\pi}}{{3{\varepsilon _{{\rm ref}}}}}$, the near-field equation becomes

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

## 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 $\textbf{g}({\textbf{r}_i},{\textbf{r}_i}) \approx - \frac{{4\pi}}{{3{\varepsilon _{{\rm 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, $\textbf{g}({\textbf{r}_i},{\textbf{r}_i}) = - \frac{{4\pi}}{{3{\varepsilon _{{\rm ref}}}}} + {\rm PV}[{\int_{{\Omega _i}} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over {\textbf{G}}}({\textbf{r}_i},\textbf{r}){\rm d}\textbf{r}}]$. Assuming the cubic subunit to be a sphere of the same volume, we get ${\rm PV}[\int_{{\Omega _i}} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over {\textbf{G}}}({\textbf{r}_i},\textbf{r}){\rm d}\textbf{r}] = \frac{2}{3}{n_{{\rm ref}}}ik_0^3$, [21]. With this better approximation, the polarizability reads [18,25]

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]

## 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]

## 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, $\textbf{E} \approx {\textbf{E}_{{\rm 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, ${\textbf{E}_{{\rm loc}}} \approx {\textbf{E}_{{\rm ref}}}$ so that $\textbf{E} \approx \frac{{3{\varepsilon _{{\rm ref}}}}}{{\varepsilon + 2{\varepsilon _{{\rm ref}}}}}{\textbf{E}_{{\rm ref}}}$. Double scattering can be introduced using the next term in the Born series, $\textbf{E} \approx {\textbf{E}_{{{{\rm Born},1}}}} = {\textbf{E}_{{\rm ref}}} + \textbf{A}{\textbf{D}_\chi}{\textbf{E}_{{\rm 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 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

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

## 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 \;{\rm 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°, and 32°. The position of the focal plane is set at ${z_o} = 2500 \; {\rm 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 _{{\rm 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 _{{\rm FFT}}}$ for multiple $(x,y)$ cuts at different $z$. The sample and illumination are the same as in 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). [CrossRef]

**2. **E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J. **186**, 705–714 (1973). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**8. **B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A **11**, 1491–1499 (1994). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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. [CrossRef]

**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). [CrossRef]

**13. **P. C. Chaumet and A. Rahmani, “Electromagnetic force and torque on magnetic and negative-index scatterers,” Opt. Express **17**, 2224–2234 (2009). [CrossRef]

**14. **N. K. Balla, P. T. So, and C. J. Sheppard, “Coupled dipole model for nonlinear scattering,” in *Frontiers in Optics 2009/Laser Science XXV/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). [CrossRef]

**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). [CrossRef]

**18. **B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. **333**, 848–872 (1988). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**21. **P. C. Chaumet, A. Sentenac, and A. Rahmani, “Coupled dipole method for scatterers with large permittivity,” Phys. Rev. E **70**, 036606 (2004). [CrossRef]

**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). [CrossRef]

**23. **E. Zubko, D. Petrov, Y. Grynko, 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). [CrossRef]

**24. **A. D. Yaghjian, “Electric dyadic Green’s functions in the source region,” Proc. IEEE **68**, 248–263 (1980). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**31. **E. Abbe, “VII.—On the estimation of aperture in the microscope,” J. R. Microsc. Soc. **1**, 388–423 (1981). [CrossRef]

**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). [CrossRef]

**33. **W. Weinstein, “The reflectivity and transmissivity of multiple thin coatings,” J. Opt. Soc. Am. **37**, 576–581 (1947). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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. [CrossRef]

**43. **R. C. Singleton, “On computing the fast Fourier transform,” Commun. ACM **10**, 647–654 (1967). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**49. **P. C. Chaumet and C. Billaudeau, “Coupled dipole method to compute optical torque: application to a micropropeller,” J. Appl. Phys. **101**, 023106 (2007). [CrossRef]

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

**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). [CrossRef]

**52. **K. Belkebir, P. C. Chaumet, and A. Sentenac, “Superresolution in total internal reflection tomography,” J. Opt. Soc. Am. A **22**, 1889–1897 (2005). [CrossRef]

**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). [CrossRef]

**54. **R. Carminati, “Phase properties of the optical near field,” Phys. Rev. E **55**, R4901–R4904 (1997). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]