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.

 

Fig. 1. (a) Object in a homogeneous background: ${\varepsilon _{{\rm ref}}} = 1$; (b) object embedded in a multilayer; the reference permittivity ${\varepsilon _{{\rm 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$.

Download Full Size | PDF

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

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.

 

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. $\textbf{k}$ is the wave vector of a plane wave diffracted by the object and $\textbf{k}^\prime $ 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.

Download Full Size | PDF

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

 

Fig. 3. Relative errors between IFDDA and Mie simulations for a nonabsorbing sphere of radius $a = 10\;\unicode{x00B5}{\rm m}$ and three different discretization steps ($N$ is the number of subunits in $\Omega_{\rm FFT}$) as a function of the sphere permittivity. The sphere is illuminated by a plane wave at $\lambda = 500 \; {\rm 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}$.

Download Full Size | PDF

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.

 

Fig. 4. Same as Fig. 3, but the sphere permittivity is complex with ${\rm Re}(\varepsilon) = 2$ and the imaginary part is increased. (a) Extinction cross section; (c) absorbing cross section; (b) optical force; and (d) asymmetry factor.

Download Full Size | PDF

 

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$, ${\rm Im}(\varepsilon) = 1$ illuminated by a plane wave as a function of ${\rm Re}(\varepsilon)$ 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}$.

Download Full Size | PDF

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.

 

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.

Download Full Size | PDF

 

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^\circ}$ with respect to the $z$ axis. (a), (c) TE polarization; (b), (d) TM polarization.

Download Full Size | PDF

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

 

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 \; {\rm 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.

Download Full Size | PDF

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.

 

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.).

Download Full Size | PDF

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

(C1)$$\textbf{E}({\textbf{r}_i}) = {\textbf{E}_{{\rm ref}}}({\textbf{r}_i}) + \sum\limits_{j = 1,i \ne j}^N \!\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over {\textbf{g}}}({\textbf{r}_i},{\textbf{r}_j}){d^3}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \chi} ({\textbf{r}_j})\textbf{E}({\textbf{r}_j})- \frac{{4\pi \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \chi} ({\textbf{r}_i})}}{{3{\varepsilon _{{\rm ref}}}}}\textbf{E}({\textbf{r}_i}).$$

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

(C2)$${\textbf{E}_{{\rm loc}}}({\textbf{r}_i}) = {\textbf{E}_{{\rm ref}}}({\textbf{r}_i}) + \sum\limits_{j = 1,i \ne j}^N \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over {\textbf{g}}}({\textbf{r}_i},{\textbf{r}_j}){\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \alpha} _{{\rm CM}}}({\textbf{r}_j}){\textbf{E}_{{\rm loc}}}({\textbf{r}_j}),$$

where

(C3)$${\textbf{E}_{{\rm loc}}}({\textbf{r}_i}) = \left({\frac{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \varepsilon} ({\textbf{r}_j}) + 2{\varepsilon _{{\rm ref}}}}}{{3{\varepsilon _{{\rm ref}}}}}} \right)\textbf{E}({\textbf{r}_i}),$$

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

(C4)$${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \alpha} _{{\rm CM}}}({\textbf{r}_j}) = \frac{3}{{4\pi}}{d^3}{\varepsilon _{{\rm ref}}}(\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \varepsilon} ({\textbf{r}_j}) - {\varepsilon _{{\rm ref}}})(\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \varepsilon} ({\textbf{r}_j}) + 2{\varepsilon _{{\rm ref}}}{)^{- 1}}$$

denotes the polarizability of each subunit related to the Clausius–Mossotti (CM) relation [18,20]. It is sometimes convenient to use this expression, as $\textbf{p} = \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over \alpha} {\textbf{E}_{{\rm 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 $\textbf{g}({\textbf{r}_i},{\textbf{r}_i})$ 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 $\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]

(D1)$${\alpha _{{\rm RR}}} = \frac{{{\alpha _{{\rm CM}}}}}{{1 - \frac{2}{3}ik_0^3{n_{{\rm ref}}}{\alpha _{{\rm CM}}}}}.$$

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]

(D2)$${\alpha _{{\rm GB}}} = \frac{{{\alpha _{{\rm CM}}}}}{{1 - \frac{2}{3}ik_0^3{n_{{\rm ref}}}{\alpha _{{\rm CM}}} - k_0^2{\alpha _{{\rm CM}}}/a}},$$

that of Lakhtakia [25]

(D3)$${\alpha _{{\rm LA}}} = \frac{{{\alpha _{{\rm CM}}}}}{{1 - 2\frac{{\varepsilon - {\varepsilon _{{\rm ref}}}}}{{\varepsilon + 2{\varepsilon _{{\rm ref}}}}}\left[{(1 - i{k_0}{n_{{\rm ref}}}a){e^{i{k_0}{n_{{\rm ref}}}a}} - 1} \right]}},$$

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

(D4)$${\alpha _{{\rm LR}}} = \frac{{{\alpha _{{\rm CM}}}}}{{1 + {\alpha _{{\rm CM}}}\left[{\frac{{({b_1} + \varepsilon {b_2}/{\varepsilon _{{\rm ref}}} + \varepsilon {b_3}/{\varepsilon _{{\rm ref}}}S)k_0^2}}{d} - \frac{2}{3}i{n_{{\rm ref}}}k_0^3} \right]}},$$

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]

(E1)$${F_u}({\textbf{r}_i}) = \frac{1}{2}\left({\sum\limits_{v = 1}^3 {p_v}({\textbf{r}_i})\frac{{\partial E_v^*({\textbf{r}_i})}}{{{\partial u}}}} \right),$$

where $u$ and $v$ stand for either $x$, $y$, or $z$ and * denotes the complex conjugate. The net force is given by $\textbf{F} = \sum\nolimits_i \textbf{F}({\textbf{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]

(E2)$$\Gamma ({\textbf{r}_i}) = {\textbf{r}^g_i} \times \textbf{F}(\textbf{r}_i) + \frac{1}{2}{\rm Re}\big[{\textbf{p}({\textbf{r}_i}) \times {{\big[{\textbf{p}({\textbf{r}_i})/{\alpha _{{\rm CM}}}({\textbf{r}_i})} \big]}^*}} \big],\!$$

where $\textbf{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\nolimits_i \Gamma ({\textbf{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.

 

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, ${\rm 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 \; {\rm nm}$; second line ${z_o} = 2500 \; {\rm nm}$; third line ${z_o} = 2800 \; {\rm nm}$.

Download Full Size | PDF

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.

 

Fig. 11. Same as Fig. 10, but the focal plane is set at ${z_o} = 2500 \; {\rm 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}$.

Download Full Size | PDF

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

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

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

(F2)$$\begin{split}\textbf{E}({\textbf{r}_\parallel},z + d) &= {e^{i{k_0}n({\textbf{r}_\parallel},z + d)d}}\\&\quad\times{\cal F}_{{2{\rm D}}}^{- 1}\left[{{{\cal F}_{{2{\rm D}}}}[\textbf{E}({\textbf{r}_\parallel},z)]({\textbf{k}_\parallel}){e^{- i({k_0} - {k_z})d}}} \right],\end{split}$$

where ${{\cal F}_{{2{\rm D}}}}$ 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).

 

Fig. 12. Modulus of the electromagnetic field on multiple $(x,y)$ cuts at different $z$ of the FFT box ${\Omega _{{\rm 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\;\unicode{x00B5}{{\rm m}^2}$, and each image has its own color scale.

Download Full Size | PDF

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]