Extended discrete dipole approximation and its application to bianisotropic media

R. Alcaraz de la Osa; P. Albella; J. M. Saiz; F. González; F. Moreno

doi:10.1364/OE.18.023865

1. Introduction

Recent advances in nanotechnology and nanoscience have been performed, involving materials with magnetic properties (magneto-optical materials) [1, 2, 3, 4, 5] as well as those with unconventional optical properties (metamaterials) [6, 7] and awakening a growing interest in this matter. The possibility of modeling the properties of such materials, and even design them in a customized way, is highly appreciated and currently constitutes a hot point. Furthermore, comparison and further agreement between experimental results and some theoretical framework is a constant requirement by researchers in those fields. In particular, the presence of a relative magnetic permeability μ_r different from 1 (even negative or tensorial) requires the revision and upgrading of the commonly available electromagnetic numerical methods [8, 9]. The widely accepted behaviour for all materials in the optical range assumes no magnetic response at high frequencies (μ_r = 1). However, effective values for the magnetic permeability different from 1 have been found for real materials within the optical range, mainly as an effect of electric currents localized within the sub-micron structure of the material, giving rise to an inherent anisotropy [10, 11]. For this reason, implementing μ_r ≠ 1 in a realistic way means considering both ɛ_r and μ_r as tensorial constituent constants. In this sense, our Group has been working, during the last years, in the generalization of some widely used numerical methods, such as the extinction theorem [12]. Other approaches include Mie theory applied to materials with arbitrary optical constants ɛ_r and μ_r [13, 14] and the discrete dipole approximation (DDA by its acronym in English) [15, 16]. In this work we focus on the second: an extension of the DDA to the case of bianisotropic materials (both ɛ_r and μ_r tensorial magnitudes), that we shall refer to as E-DDA, and is described in the next section. Commercial software, like COMSOL and FDTD, can also deal with bianisotropic media. In fact, they have been used in many researches to solve 1-D [19] and 2-D [20] problems. In this work, where 3-D systems are involved, DDSCAT (based on DDA) [17, 18], has been used, when possible, to check the reliability and precision of our calculations, being very helpful in the understanding of the required new point of view.

2. E-DDA code

2.1. Method description

E-DDA stands for Extended Discrete Dipole Approximation. The formalism of the currently available DDA methods is restricted to materials with relative magnetic permeability μ_r = 1. One of the main contributions of this work is the generalization of the method to materials with arbitrary electric and magnetic susceptibility tensors (bianisotropic materials). Along this section, we will describe in detail the theoretical formalism of the E-DDA method. As You et. al [8] have already suggested, following Lakhtakia’s steps, at each lattice site we must locate two dipoles, one electric and one magnetic, so that we can take into account both the electric and magnetic responses of the material. The dipole moment p of an electric dipole under the action of an electric field E is given by p = ɛ₀ɛ_mα̿E, where ɛ_m is the relative electric permittivity of the surrounding media, and α̿ is the electric polarizability tensor. In the same way, the dipole moment m of a magnetic dipole under the influence of a magnetic field H is given by $m = \frac{\bar{\bar{χ}}}{μ_{0} μ_{m}} B = \bar{\bar{χ}} H$ , where μ_m is the relative magnetic permeability of the surrounding media, and χ̿ is the magnetic susceptibility tensor. Both α̿ and χ̿ can be related to the optical properties (ɛ̿_r, μ̿_r) through the well-known Clausius-Mossotti relation:

{\bar{\bar{α}}}_{CM} = 3 V ({\bar{\bar{ɛ}}}_{r} - ɛ_{m} \bar{\bar{I}}) {({\bar{\bar{ɛ}}}_{r} + 2 ɛ_{m} \bar{\bar{I}})}^{- 1}

{\bar{\bar{χ}}}_{CM} = 3 V ({\bar{\bar{μ}}}_{r} - μ_{m} \bar{\bar{I}}) {({\bar{\bar{μ}}}_{r} + 2 μ_{m} \bar{\bar{I}})}^{- 1}

A radiative correction is also included to take into account the phase lag between the incident light and the scattered light radiated by the dipole [21]:

\bar{\bar{α}} = {\bar{\bar{α}}}_{CM} {(\bar{\bar{I}} - \frac{i k^{3} {\bar{\bar{α}}}_{CM}}{6 π})}^{- 1}

\bar{\bar{χ}} = {\bar{\bar{χ}}}_{CM} {(\bar{\bar{I}} - \frac{i k^{3} {\bar{\bar{χ}}}_{CM}}{6 π})}^{- 1}

Although this radiative correction is precise enough in most of the real situations, further corrections on the dipole polarizabilities are being used by several authors, developing prescriptions aimed to improve the accuracy [22, 23].

As in the conventional DDA [17, 24], each dipole feels an electromagnetic field sum of the incident field at that site, and the contributions due to all the other dipoles. Working with the same formalism used by Mulholland et. al [25], which utilizes the M.K.S. unit’s system, we can express the total electric and magnetic fields at lattice site j as:

E_{j} = E_{inc, j} + \sum_{k \neq j}^{N} {\bar{\bar{α}}}_{j} {\bar{\bar{C}}}_{j k} E_{k} - \sqrt{\frac{μ_{0} μ_{m}}{ɛ_{0} ɛ_{m}}} \sum_{k \neq j}^{N} {\bar{\bar{χ}}}_{j} {\bar{\bar{f}}}_{j k} H_{k}

H_{j} = H_{inc, j} + Σ_{k \neq j}^{N} {\bar{\bar{χ}}}_{j} {\bar{\bar{C}}}_{j k} H_{k} + \sqrt{\frac{ɛ_{0} ɛ_{m}}{μ_{0} μ_{m}}} Σ_{k \neq j}^{N} {\bar{\bar{α}}}_{j} {\bar{\bar{f}}}_{j k} E_{k}

, where E_inc,_j = E₀ exp(ik · r_jk) and H_inc,_j = H₀ exp(ik · r_jk), with

H_{0} = {[ɛ_{0} ɛ_{m} / (μ_{0} μ_{m})]}^{1 / 2} (k \times E_{0}^{*})

and where label[*] denotes the complex conjugate. The matrixes C̿_jk and f̿_jk involved in Eqs. (5) and (6) can be found in Ref. [25].

In order to solve the set of equations in Eqs. (5) and (6), it is convenient to rewrite them using its associated matrix form. Then, one can easily obtain the system of equations in the form:

\bar{\bar{A}} x = b

, where A̿ is a 6N × 6N matrix, x is the 6N-dimensional vector of unknowns given by x = (E₁, H₁, E₂, H₂, . . . , E_N, H_N) and b is the 6N-dimensional vector of independent terms containing the incident field at each lattice site: b = (E_inc,1, H_inc,1, E_inc,2, H_inc,2, . . . , E_inc,_N, H_inc,_N). This system can be solved using some of the iterative methods already available. Initially we used the successive approximations method, proving to be very slow and not convergent on most cases. At this time, E-DDA implements the Complex-Conjugate Gradient (CCG) method through the Botchev’s subroutine [26].

2.2. Code Performance

The main achievement of our code in its current version is versatility, in the sense that it is able to deal with situations involving arbitrary electric and magnetic susceptibilities in a broad sense. In the following section we shall survey the most representative situations at reach for this code, from the most basic one of a dielectric, to the most generals, like bianisotropic materials, metals, metamaterials (in particular left-handed materials, with both real parts of ɛ_r and μ_r negative) or magneto-optical materials (with ɛ̿_r an antisymmetric tensor). Some of these situations admit direct comparison with past DDA versions, while other can only be treated with other calculations methods, and not in a feasible way. The former will contribute to validate our code, while the latter constitute genuine novel results. Both together, prove the potential of E-DDA as a new and reliable computing tool. It is interesting to remark that the incident wave polarization can be arbitrary (elliptical in the most general case), allowing polarimetric calculations to be performed. As for extinction and absorption cross-sections (and efficiencies), they have been implemented following [9]:

C_{ext} = \frac{k}{ɛ_{0} ɛ_{m} {| E_{0} |}^{2}} \sum_{j = 1}^{N} ℑ [E_{inc, j}^{*} \cdot p_{j} + μ_{0} μ_{m} H_{inc, j}^{*} \cdot m_{j}]

C_{abs} = \frac{k}{ɛ_{0} ɛ_{m} {| E_{0} |}^{2}} \sum_{j = 1}^{N} {ℑ [E_{j}^{*} \cdot p_{j}] - \frac{k^{3}}{6 π ɛ_{0} ɛ_{m}} {| p_{j} |}^{2} + μ_{0} μ_{m} [ℑ [H_{j}^{*} \cdot m_{j}] - \frac{k^{3}}{6 π} {| m_{j} |}^{2}]}

The scattering cross-section can be easily obtained by the difference of the extinction and absorption cross-sections: C_sca = C_ext – C_abs, but a more convenient way is to compute the far-field scattered by the object [27]:

C_{sca} = \frac{k^{2}}{16 π^{2} ɛ_{0} ɛ_{m} {| E_{0} |}^{2}} {\int | \sum_{j = 1}^{N} exp (- i k n \cdot r_{j}) {\frac{1}{\sqrt{ɛ_{0} ɛ_{m}}} [p_{j} - (n \cdot p_{j}) n] - \sqrt{μ_{0} μ_{m}} n \times m_{j}} |}^{2} d Ω

, where n is an unit vector in the direction of scattering. Additionally, we can also define the phase-lag cross-section in terms of the imaginary part of the forward-scattering amplitude [21]:

C_{pha} = \frac{k}{2 ɛ_{0} ɛ_{m} {| E_{0} |}^{2}} \sum_{j = 1}^{N} ℜ [E_{inc, j}^{*} \cdot p_{j} + μ_{0} μ_{m} H_{inc, j}^{*} \cdot m_{j}]

Conventional angular scattering patterns (scattered intensities) can be obtained, as well as the full Mueller Matrix of the system, completing the far-field results list. As for local magnitudes, these include dipole moments, fields or Poynting vector distributions, local phase functions, etc.

3. Testing the code

To check the reliability of E-DDA, we have performed some calculations on systems that conventional methods can solve. Fig. 1 shows the extinction, absorption and scattering efficiencies for both a gold sphere of radius R = 20 nm, and a sphere with a dielectric (ɛ_c = 2) core (inclusion) of radius R = 12 nm and a metallic (gold, optical constants taken from Johnson and Christy [28]) shell, for an external radius R = 20 nm. Comparison between our code and the well-proved DDSCAT code from B. T. Draine [17, 18] is also presented, finding a very good agreement in both, the spectral shape and the absolute differences for all the efficiencies within the optical range.

Fig. 1 Extinction, absorption and scattering efficiencies for both a gold sphere of radius R = 20 nm, and a sphere with a dielectric (ɛ_c = 2) core (inclusion) of radius R = 12 nm and a metallic (gold) shell, for an external radius R = 20 nm. Comparison between our code and the well-established DDSCAT code is also provided. The dipole spacing was d = 4 nm, with N = 515.

Download Full Size | PDF

4. Bianisotropic media

When the E-DDA is used in the most general case of a bianisotropic scatterer, it is possible to find a situation in which we already have a knowledge of the system’s behavior. Such is the case of particles with directional scattering imposed by their optical constants as, for instance, the null-scattering conditions proposed by Kerker et. al. [29]. It will be shown here that these well-known conditions remain valid for bianisotropic media. As a reminder, Kerker’s null scattering conditions are obtained when the electric and magnetic polarizabilities meet certain stipulations. In particular, when α̿ = χ̿ the backscatter gain equals zero, and when α̿ = –χ̿ the forward scatter is zero. These conditions lead to ɛ̿_r = μ̿_r for the first case, and ɛ̿_r = (4I̿ – μ̿_r)(2μ̿_r + I̿)⁻¹ for the latter. However, these relations fail to include the radiative correction [21]. It can be shown that the zero-backward scattering condition remains valid, that is, ɛ̿_r = μ̿_r. But, for the zero-forward scattering condition, we obtain the new relation:

{\bar{\bar{μ}}}_{r} = [4 \bar{\bar{I}} - {\bar{\bar{ɛ}}}_{r} - \frac{i {(k d)}^{3}}{π} ({\bar{\bar{ɛ}}}_{r} - \bar{\bar{I}})] {[2 \bar{\bar{I}} + {\bar{\bar{ɛ}}}_{r} - \frac{i {(k d)}^{3}}{π} ({\bar{\bar{ɛ}}}_{r} - \bar{\bar{I}})]}^{- 1}

Let us now consider a set of three different permittivity tensors. The first case (isotropic material) is the scalar case:

{\bar{\bar{ɛ}}}_{r_{1}} = (\begin{matrix} 2.0 + 0.01 i & 0 & 0 \\ 0 & 2.0 + 0.01 i & 0 \\ 0 & 0 & 2.0 + 0.01 i \end{matrix}) = (2.0 + 0.01 i) \bar{\bar{I}}

The second case is chosen in the form of a typical magneto-optical material, with an antisymmetric relative electric permittivity tensor ɛ̿_r₂ given by (magnetization along z-axis):

{\bar{\bar{ɛ}}}_{r_{2}} = (\begin{matrix} 2.0 + 0.01 i & 0.3 + 0.2 i & 0 \\ - 0.3 - 0.2 i & 2.0 + 0.01 i & 0 \\ 0 & 0 & 2.0 + 0.01 i \end{matrix})

For the third case we propose a symmetric permittivity tensor ɛ̿_r₃:

{\bar{\bar{ɛ}}}_{r_{3}} = (\begin{matrix} 2.0 + 0.01 i & 0.3 + 0.2 i & 0 \\ 0.3 - 0.2 i & 2.0 + 0.01 i & 0 \\ 0 & 0 & 2.0 + 0.01 i \end{matrix})

For each of these three cases we shall define two values of the relative magnetic permeability tensor in order to fulfill each of the zero scattering conditions (this makes six different materials as a whole). We now consider an sphere of diameter D = 20 nm made of those materials, illuminated with a wavelength of λ = 500 nm, obtaining that the zero-forward scattering condition is, for each case:

{\bar{\bar{μ}}}_{r_{1}} = (0.4 - 3.6005 \times 10^{- 3} i) \bar{\bar{I}}

{\bar{\bar{μ}}}_{r_{2}} = (\begin{matrix} 3.9232 \times 10^{- 1} - 2.0508 \times 10^{- 2} i & - 1.0899 \times 10^{- 1} - 6.8489 \times 10^{- 2} i & 0 \\ 1.0899 \times 10^{- 1} + 6.8489 \times 10^{- 2} i & 3.9232 \times 10^{- 1} - 2.0508 \times 10^{- 2} i & 0 \\ 0 & 0 & 0.4 - 3.6005 \times 10^{- 3} i \end{matrix})

{\bar{\bar{μ}}}_{r_{3}} = (\begin{matrix} 4.0712 \times 10^{- 1} + 1.3870 \times 10^{- 2} i & - 1.0804 \times 10^{- 1} - 7.3802 \times 10^{- 2} i & 0 \\ - 1.0804 \times 10^{- 1} - 7.3802 \times 10^{- 2} i & 4.0712 \times 10^{- 1} + 1.3870 \times 10^{- 2} i & 0 \\ 0 & 0 & 0.4 - 3.6005 \times 10^{- 3} i \end{matrix})

and μ̿_{r_i} = ɛ̿_{r_i} (i = 1, 2, 3) for the zero backward. The scattering patterns show that Kerker’s conditions are being satisfied (Fig. 2). It is worth noticing that, while the isotropic case can be computed by means of conventional numerical methods (producing a perfect match), the bianisotropic ones do not admit such comparison in a feasible way. By using E-DDA, we obtain results that show that the systems behave exactly as expected from the theoretically imposed condition, that is, exhibiting zero-back and zero-forward scattering.

Fig. 2 a) Zero-backward scattering (μ̿_r = ɛ̿_r). b) Zero-forward scattering. In every case, the dipole spacing was d = 2 nm, with N = 515.

Download Full Size | PDF

5. Conclusions

We have extended the applicability range of the discrete dipole approximation to the case of anisotropic and magnetic materials, including the bianisotropic case. We have tested the validity of the method for a case where its applicability range overlaps with the one of a well-proved code, which includes inhomogeneities and presence of metallic and dielectric media, finding a very good agreement. We have applied the proposed method to a situation out of reach for current implementations of the DDA, like the null-scattering (backward and forward) conditions for bianisotropic media. We have verified numerically these conditions in its tensorial form.

Acknowledgments

This research has been supported by MICINN (Ministerio de Ciencia e Innovacion) under project #FIS2007-60158. The authors thankfully acknowledge the computer resources provided by the RES (Red Espaola de Supercomputacion) node (Altamira) at Universidad de Cantabria. R. Alcaraz de la Osa also thanks the Ministry of Education of Spain for his FPU grant.

References and links

1. S. Albaladejo, R. Gómez-Medina, L. S. Froufe-Pérez, H. Marinchio, R. Carminati, J. F. Torrado, G. Armelles, A. García-Martín, and J. J. Sáenz, “Radiative corrections to the polarizability tensor of an electrically small anisotropic dielectric particle,” Opt. Express 18(4), 3556–3567 (2010). [CrossRef] [PubMed]

2. B. Sepúlveda, J. B. González-Díaz, A. García-Martín, L. M. Lechuga, and G. Armelles, “Plasmon-induced magneto-optical activity in nanosized gold disks,” Phys. Rev. Lett. 104(14), 147401 (2010). [CrossRef] [PubMed]

3. G. Ctistis, E. Papaioannou, P. Patoka, J. Gutek, P. Fumagalli, and M. Giersig, “Optical and magnetic properties of hexagonal arrays of subwavelength holes in optically thin cobalt films,” Nano Lett. 9(1), 1–6 (2008). [CrossRef] [PubMed]

4. D. A. Smith and K. L. Stokes, “Discrete dipole approximation for magneto-optical scattering calculations,” Opt. Express 14(12), 5746–5754 (2006). [CrossRef] [PubMed]

5. N. B. Piller and O. J. F. Martin, “Extension of the generalized multipole technique to three-dimensional anisotropic scatterers,” Opt. Lett. 23(8), 579–581 (1998). [CrossRef]

6. V. Agranovich and Y. Gartstein, “Electrodynamics of metamaterials and the Landau–Lifshitz approach to the magnetic permeability,” Metamaterials 3(1), 1–9 (2009). [CrossRef]

7. A. Alu and N. Engheta, “The quest for magnetic plasmons at optical frequencies,” Opt. Express 17(7), 5723–5730 (2009). [CrossRef] [PubMed]

8. Y. You, G. W. Kattawar, P.-W. Zhai, and P. Yang, “Zero-backscatter cloak for aspherical particles using a generalized DDA formalism,” Opt. Express 16(3), 2068–2079 (2008). [CrossRef] [PubMed]

9. P. C. Chaumet and A. Rahmani, “Coupled-dipole method for magnetic and negative-refraction materials,” J. Quant. Spectrosc. Radiat. Transf. 110(1–2), 22–29 (2009).

10. A. Alù, A. Salandrino, and N. Engheta, “Negative effective permeability and left-handed materials at optical frequencies,” Opt. Express 14(4), 1557–1567 (2006). [CrossRef] [PubMed]

11. R. Marqués, F. Medina, and R. Rafii-El-Idrissi, “Role of bianisotropy in negative permeability and left-handed metamaterials,” Phys. Rev. B 65(14), 144440 (2002). [CrossRef]

12. P. Albella, F. Moreno, J. M. Saiz, and F. González, “Backscattering of metallic microstructures with small defects located on flat substrates,” Opt. Express 15(11), 6857–6867 (2007). [CrossRef] [PubMed]

13. B. García-Cámara, F. Moreno, F. González, J. M. Saiz, and G. Videen, “Light scattering resonances in small particles with electric and magnetic properties,” J. Opt. Soc. Am. A 25(2), 327–334 (2008). [CrossRef]

14. B. García-Cámara, F. González, F. Moreno, and J. M. Saiz, “Exception for the zero-forward-scattering theory,” J. Opt. Soc. Am. A 25(11), 2875–2878 (2008). [CrossRef]

15. P. Albella, F. Moreno, J. M. Saiz, and F. González, “Surface inspection by monitoring spectral shifts of localized plasmon resonances,” Opt. Express 16(17), 12,872–12,879 (2008). [CrossRef]

16. P. Albella, J. M. Saiz, J. M. Sanz, F. González, and F. Moreno, “Nanoscopic surface inspection by analyzing the linear polarization degree of the scattered light,” Opt. Lett. 34(12), 1906–1908 (2009). [CrossRef] [PubMed]

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

18. B. T. Draine and P. J. Flatau, User Guide for the Discrete Dipole Approximation Code DDSCAT 7.1 (2010). http://arxiv.org/abs/1002.1505.

19. A. Akyurtlu and D. Werner, “Modeling of transverse propagation through a uniaxial bianisotropic medium using the finite-difference time-domain technique,” IEEE Trans. Antennas Propag. 52, 3273–3279 (2004). [CrossRef]

20. Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljacić, “Observation of unidirectional backscattering-immune topological electromagnetic states,” Nature 461(7265), 772–775 (2009). [CrossRef] [PubMed]

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

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

23. D. Gutkowicz-Krusin and B. T. Draine,“Propagation of electromagnetic waves on a rectangular lattice of polarizable points,” (2004). arXiv:astro-ph/0403082v1.

24. M. Yurkin, V. Maltsev, and A. Hoekstra, “The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength,” J. Quant. Spectrosc. Radiat. Transf. 106(1–3), 546–557 (2007). [CrossRef]

25. G. W. Mulholland, C. F. Bohren, and K. A. Fuller, “Light scattering by agglomerates: coupled electric and magnetic dipole method,” Langmuir 10(8), 2533–2546 (1994). [CrossRef]

26. M. A. Botchev, SUBROUTINE ZBCG2, http://www.math.uu.nl/people/vorst/zbcg2.f90 (2001).

27. J. D. Jackson, Classical Electrodynamics Third Edition, 3rd ed. (Wiley, 1998).

28. P. Johnson and R. Christy, “Optical constants of the noble metals,” Phys. Rev. B. 6, 4370–4379 (1972). [CrossRef]

29. M. Kerker, D.-S. Wang, and C. L. Giles, “Electromagnetic scattering by magnetic spheres,” J. Opt. Soc. Am. 73(6), 765–767 (1983). [CrossRef]

Extended discrete dipole approximation and its application to bianisotropic media

Abstract

1. Introduction

2. E-DDA code

2.1. Method description

2.2. Code Performance

3. Testing the code

4. Bianisotropic media

5. Conclusions

Acknowledgments

References and links

Cited By

Figures (2)

Equations (18)

Optics Express