The Discrete Dipole Approximation (DDA) formalism has been generalized to materials with permeabilities µ ≠ 1 to study the scattering properties of impedance-matched aspherical particles and cloaked spheres. We have shown analytically that any impedance-matched particle with a four-fold rotational symmetry with respect to the direction of the incident radiation has the feature of zero backscatter. Moreover, an impedance-matched coat with the aforementioned symmetry property acting on an irregular dielectric particle with the same symmetry property can substantially reduce the backscatter. This leads to a substantial reduction of the signals from an object being detected by a monostatic radar/lidar system. The DDA simulation also provides accurate information about electric field distributions in the vicinity of a cloaked sphere.
©2008 Optical Society of America
Invisibility cloaking based on metamaterials [1, 2, 3, 4, 5] has been an active research topic recently. Based on coordinate transformations, it has been demonstrated by Pendry et al.  that an anisotropic and inhomogeneous spherical shell gives an absolute zero scattering cross section if the permittivity and permeability tensors are rotationally uniaxial and vary as follows
are the tangential components, and
are the radial components of the permittivity and permeability, respectively, with ε0 and µ0 the free space permittivity and permeability, respectively, and R1 and R2 the inner and outer radii of the spherical cloak, respectively. Ray-tracing simulations [1, 2] supporting the conclusions by Pendry et al.  were reported in the geometric optics limit. Full-wave finite-element simulations  were performed to study the effects of cloaking material perturbations to the propagation of the incidentwaves in a 2-D cylindrical case. A 2-D cloaking has been practically realized  using artificially structured metamaterials at microwave frequencies. A rigorous solution to Maxwell’s equations for a cloaked sphere has been reported  confirming the above result.
An invisibility cloak using a homogeneous coating has also been investigated. In the context of the Lorenz-Mie theory, Alù et al.  explored the possibilities of substantial reduction of the scattering cross section of small spherical objects and cylindrical objects using homogeneous coats with negative or low permittivity/permeability. Zhai et al.  noticed that the cloaking condition is not very strict if only backscatter is concerned. They showed that when the impedance was equal to the impedance of the surrounding medium, Z = Z0, in both inner and outer spheres with homogeneous materials, an exact zero backscatter is guaranteed in the case of a coated sphere. Note that an analog to this phenomenon, impedance matching, is well known to electrical engineers. When any form of energy is transferred between a source and a load, making the output impedance of the source equal to the input impedance of the load maximizes the power transfer and minimizes reflections from the load. It also applies to the transfer of light between two dielectric media. This concept, however, has drawn little attention in the light scattering community. Here, to obtain a zero backscatter, the impedance of the scattering particle is set to match the impedance of the surrounding medium. The latter is obviously unity. The authors of Ref. further showed that, even if the inner core does not satisfy the impedance-matching condition, the backscatter could still be very small, provided that the outer coat is impedance-matched and the product of the thickness and the imaginary part of the refractive index of the outer coat is sufficiently large. This observation has an important implication to monostatic radar/lidar detection.
To date, analytical solutions are available for regular-shaped particles only. On the other hand, numerical approaches, such as the ray-tracing technique and the finite-element method, can be readily applied to aspherical particles.
The Discrete Dipole Approximation (DDA) [8, 9, 10] method is a robust numerical technique for computation of elastic scattering and absorption by dielectric particles with arbitrary shapes. In the DDA, a dielectric particle is represented by a cubic array of electric dipoles, each of which oscillates in response to the incident radiation and the scattered field induced by all the other dipoles. The Clausius-Mossotti relation  is normally used to relate the polarizabilities of these dipoles to the local dielectric constant. Consequently, a self-consistent solution for all the dipole moments can be obtained, and thus the scattered electric field can be determined.
In this study, we will first present a generalization of the conventional DDA formalism applicable to materials with permeabilities µ ≠ 1. Next we will use this formalism to show analytically that any impedance-matched particle with four-fold rotational symmetry has zero backscatter. Then we will compare the results obtained from the DDA simulations with the analytical results for coated spheres. We will further confirm the above statement on zero backscatter by numerical simulations. Finally, we will present simulated near- and far-fields in the case of cloaked spheres.
2. Generalization of the DDA method
Applications of the conventional DDA formalism are limited to dielectric particles with permeability µ = 1 and vanishing magnetic susceptibility χm=µ − 1 = 0, as the null magnetization M = χmH = 0 allows us to represent the particle solely by a collection of electric dipoles. For materials with nonzero magnetic susceptibilities, e.g., bianisotropic materials and metamaterials, the magnetization M is non-vanishing. In this case, as has been reported by Lakhtakia , a collection of both electric and magnetic dipoles should be used to account for both electric and magnetic responses of the particle to the incident electromagnetic radiation.
Following Lakhtakia’s steps, we locate an electric and a magnetic dipole at each lattice site j in the DDA cubic lattice with j = 1, …, N running over all occupied lattices sites. Each electric dipole is characterized by a polarizability tensor α̿j such that Pj=α̿j E ext,j, where Pj is the instantaneous electric dipole moment and E ext, j is the instantaneous electric field at position j due to all sources external to j, i.e., the incident radiation and the oscillating electric and magnetic dipoles at all other N-1 locations. Similarly, each magnetic dipole is characterized by a magnetic susceptibility tensor χ̿ m, j such that M j=χ̿mj Hext, j, where Mj is the instantaneous magnetic dipole moment and Hext, j is the instantaneous magnetic field at position j due to all sources external to j.
Like the Clausius-Mossotti equation, for paramagnetic materials there is also an equation  that relates the microscopic magnetic susceptibility for each magnetic dipole χ̿ m, j and the local macroscopic permeability µ̿j in the zero frequency limit kd→0, where k is the wavenumber of the incident raidation, and d is the inter-dipole spacing. In the coordinate system where the permeability and susceptibility tensors are diagonalized, the relation is
for each diagonal component ll, where n=N/V is the number density of dipoles, with N the number of dipoles and V the volume of the particle. In the long-wavelength limit, kd≪1, the lattice dispersion relation (LDR)  also holds for the magnetic susceptibility, giving a magnetic susceptibility at finite frequency χ̿m(kd).
Knowing the polarizability α̿j and the magnetic susceptibility χ̿ m, j at all lattice sites j, the self-consistent set of equations for electric and magnetic dipole moments P j and M j (j=1, …, N) can be written as
where in the first term
are the electric and magnetic fields at position j due to the incident electromagnetic radiation, with E0 × H0 = k, the wave vector of the incident radiation. In the second term, and are the contributions to the electric field at position j due to the electric and magnetic dipoles at position k, given by [9, 11]
and and are the contributions to the magnetic field at position j due to the electric and magnetic dipoles at position k, given by  (pp. 413)
where rjk = rj − r k and rjk = |rjk|. Notice here that the electric dipole moments and magnetic dipole moments are coupled with each other.
Let Ajk indicate the off-diagonal elements of A-matrices, and define the diagonal elements of these A-matrices as follows
Then, Eq.(5) can be reformulated into a compact form as a set of 2N inhomogeneous linear complex vector equations
Note here that each element of the A-matrices, Ajk, is a 3×3 matrix and that and are symmetric, while and are anti-symmetric with respect to lattice site indices i and j. Therefore, we can still use the Fast Fourier Transform (FFT) and the Complex-Conjugate Gradient (CCG) methods used in the conventional DDA formalism to solve for P j’s and Mj’s in Eq.(10). The CCG algorithm takes more iterations to converge here, since a 6N×6N matrix instead of a 3N×3N matrix is involved.
3. Backscatter of impedance-matched particles with four-fold rotational symmetry
The above equations describe the responses of a particle to the incident electric and magnetic radiation in a discretized fashion. They can be readily used to study the backscatter of impedance-matched particles with four-fold rotational symmetry. Here, we limit our discussions to situations where the particle is embedded in a vacuum, in which case the impedancematching condition becomes ε = µ for the particle. Let us consider an arbitrary particle with four-fold rotational symmetry with respect to the incident radiation, as shown in Fig. 1(I). We assume that the incident radiation is along the z-axis.
The local electric and magnetic moments at 1 can be denoted with the following components
Now we rotate the coordinate system clockwise by 90°. The particle positions look as shown in Fig. 1(II) in the rotated coordinate system. This is equivalent to rotating the particle and the fields counterclockwise with the coordinate system fixed. Since this operation is only a rotation of the coordinate system, the electromagnetic field remains the same. Therefore, the local electric and magnetic moments at 1 will not change. In the rotated coordinate system, they read
P1′ = (−y,x,z), M1′ = (−ν,u,w).
Next, we relabel the electric and magnetic fields of the incident radiation as follows
E″ = −H′, H″ = E′.
Here, we do not attach any physical meaning to the quantities E″ and H″. This operation turns Eq.(10) into
or in a rearranged form
For a particle with unit impedance α̿j=χ̿ m, j, the following relationships hold for any pair of indices j and k:
Thus, Eq.(12) can be rewritten as
A comparison of Eq.(13) with Eq.(10) immediately reveals that E″ and H″ can be interpreted as electric and magnetic parts of an incident radiation, and the local electric and magnetic moments as responses to E″ and H″ can be readily recognized as follows:
which implies that
P 1″ = (ν,−u,−w), M 1″ = (−y,x,z).
Following the same procedure, we can get
The scattered electric field due to the electric and magnetic dipoles at these four points can then be written as 
where n is the unit vector along a scattering direction. At the backscatter direction where n = (0,0,−1), it is straight forward to check that E1234(180°) vanishes. For any particle with four-fold rotational symmetry with respect to the incident radiation, the total scattered electric field at 180° can be written as a superposition of a collection of E1234’s, therefore the backscatter is always zero for such a particle.
This interesting result attributes to a simple fact that the electric and magnetic parts of the incident radiation are perpendicular to each other, and that the particle has equal electric permittivity and magnetic permeability. Therefore, the electric and magnetic responses of the particle can be switched under four-fold rotations.
4. Numerical results
To validate the accuracy of the generalized DDA formalism, we compared the DDA results with the analytical counterpart obtained by the Lorenz-Mie theory  for coated spheres with matched impedances. Here for both layers the refractive index is , and the impedance is . In our comparison, we used the following parameters: the inner core is a regular dielectric medium with n1 = 1.5 and Z1 = 1/1.5, or and ; the outer coating is an absorbing medium with unit impedance, n2 = 2+i and Z2 = 1, or , with being the imaginary unit. The size parameters for the inner and outer layers are denoted by x=2πR 1/λ and y=2πR 2/λ, respectively. An inter-dipole spacing d~λ/40 was used in the DDA simulations. The resulting P11, the (1,1) element of the phase matrix, computed from the Lorenz-Mie theory and the DDA formalism are shown in Fig. 2. In spherical cases, P11 is identical in scattering planes with any azimuthal angle ϕ due to symmetry. This figure shows that the DDA results agree with the Lorenz-Mie results very well for coated spheres of size parameter as large as y = 8.
We then explored the scattering of a cube with unit impedance. Here we characterize the size of an aspherical particle in terms of the size parameter of an equal-volume sphere x = 2πa eff/λ where aeff = (3V/4π)1/3 with V the volume of the aspherical particle. Shown in Fig. 3 is P11 for a cube of size parameter x = 8 in scattering planes ϕ = 0° and ϕ = 45°. The cube is oriented with one face facing the incident radiation. Thus, there is a four-fold rotational symmetry with respect to the direction of the incident radiation. The scatterer has a refractive index of n = 1.5 and an impedance of Z = 1, or . For comparison, results for a regular dielectric particle (n = 1.5 and Z = 1/1.5, or and ) with the same shape are also shown.
Figure 3 first implies that the size parameter x = 8 is sufficient to show effects of asphericity, as P11 in the two scattering planes are appreciably different. It turned out that for an impedance-matched cubic particle in this specific orientation, the simulated P11 remains almost the same as that for a regular dielectric particle in forward scattering directions, and starts to differ from the dielectric particle result at scattering angle θ > 90°. It starts to drop at about θ = 170°, and is on the order of 10-9 at 180°, the same order as the simulated backscatter for spheres. It is reasonable to say that this very small non-vanishing value at 180° is due to numerical errors, and the backscatter in this case remains exactly zero as predicted before.
We also did simulations for the following columns with one base facing the incident radiation: an octagonal column with an eight-fold rotational symmetry; a hexagonal column with a six-fold symmetry; and a cuboid, or a rectangular column, with a two-fold symmetry. We kept the volume the same as before with a size parameter x = 8, and the aspect ratio (the length of the column to the diameter of the base of the column h/d) equal to one. The projections of these columns on the x-y plane are shown in Fig. 4(a). In our simulations, we set n = 1.21 and Z = 1 to represent a soft impedance-matched medium. The resulting P11’s are shown in Fig. 4(b), along with the results for a cube.
Figure 4 shows that, other than the cube, the octagonal column is the only one that has the zero backscatter feature. Noticing that the octagonal column also has a four-fold rotational symmetry, this result is not surprising. These numerical simulations confirm our previous prediction that any impedance-matched particle with a four-fold rotational symmetry with respect to the incident radiation has zero backscatter.
Moreover, we show the resultant scattering for a coated cube, where a smaller cube is embedded in a larger cubic shell with their centers at the same location. Similarly to the coated sphere case, the regular dielectric inner cube has a size parameter x = 4 and a non-unit impedance (n1 = 1.1 and 1.5, Z1 = 1/n1), and the absorbing outer cubic coating has a size parameter y=8 and a unit impedance (n2 = 2+i and Z2 = 1). The results are shown in Fig. 5, where the blue curve corresponds to the n1 = 1.1 case, and the red curve corresponds to the n1 = 1.5 case. The dash curves show the scattering of light by uncoated dielectric cores for comparison.
As can be seen in Fig. 5, the backscatter of a coated cube is substantially reduced compared with that of the corresponding uncoated cubic core. The resulting backscatter is as low as ~10−7, which is of the same order of the simulated backscatter in the coated-sphere case shown above. What is more, other than a minor difference at the backscatter direction, the whole scattering pattern is almost insensitive to what is located in the coated region.
The DDA formalism can also be applied to the perfect cloak with tensorial permittivity and permeability . The difficulty for DDA, however, is that in the cloak material, the tangential component of the refractive index determined by Eq.(2) could be very large if the ratio R 1/R2 is close to 1. For example, R1 = 0.3R2 requires nt = 1.43, which is a moderate value; while R1 = 0.5R2 requires nt = 2, which could be prohibitive for the DDA formalism. Moreover, the inhomogeneous cloak represented by the discrete array of dipoles is not a perfect cloak. Therefore, the cloaking effects are not expected to be perfect.
Shown in Fig. 6 is the electric field in the vicinity of a cloaked sphere with a cloak described by Eq.(1). In the simulations, we chose R2 = 8λ/2π, and the results for two cases, R1 = 0.3R2 and 0.5R2 are shown. For comparison, we also show the corresponding field distribution associated with an uncloaked homogeneous impedance-matched sphere with n = 1.21 and Z = 1. A unit amplitude was assumed for the incident electric field.
A quick comparison between the cloaked and uncloaked results reveals that the DDA calculation simulates the plane-wave fields outside of the cloaked particle and the distorted fields in the cloaking material correctly. The DDA simulation of the electric field inside the cloaked region, however, is not exactly zero. In the R1 = 0.3R2 case, a maximum leakage of about 13% of the radiating field, or about 1.7% of the radiative energy, into the cloaked region is observed due to the discretization of material properties in the DDA simulation. In the R1=0.5R2 case, the maximum leakage is about 24% of the radiating field, or about 6% of the radiative energy.
The same electric field distribution associated with cloaked spheres can be found in previous studies [3, 5], where Cummer et al.  presented a field distribution simulated by the finite-element method, and Chen et al. showed the analytical results by solving Maxwell’s equations. Comparing the DDA simulation with these results, we found that the DDA gives better plane-wave structure outside of the cloak than the finite-element method does.
As a numerical method, the DDA is not expected to give an exact zero scattering cross section of the cloaked sphere. However, the simulated scattering efficiencies for the cases shown in Fig. 6 are on or below the order of 10-4 (6.13×10-5 in case (a) with R1 = 0.3R2 and 4.43×10-4 in case (b) with R1 = 0.5R2), which are sufficiently small to exhibit the nature of invisibility. For comparison, the scattering efficiency of the uncloaked sphere in case (c) is Qsca = Csca/(πR2) = 3.46. We use the outer radius R2 in the calculation of Qsca as we are comparing cases with the same R2. It can also be noticed that the DDA calculations gave a better simulation of the zero scattering for the cloak with R1 = 0.3R2, in which case nt is smaller.
To study the capability of the DDA formalism in simulating the far field scattered from a cloaked object, in Fig. 7 we show the simulated P11 for a cloaked sphere with R2 = 5λ/2π, R1 = 0.2R2, and various refractive indices in the cloaked region. The P11 for uncloaked cores are also shown for comparison. The corresponding scattering efficiencies are listed in Table 1. These results imply that the DDA formalism has its limitation in simulating the zero scattering of cloaked spheres, as the simulated P11 and Qsca of the cloaked spheres are only 3 orders lower than that of the uncloaked cores. However, as the refractive index of the cloaked region increases and the scattering of the uncloaked core increases, the simulated scattering of the cloaked spheres stays almost the same. This is a demonstration of the cloaking effects.
The Discrete Dipole Approximation method is generalized to materials with permeabilities µ ≠ 1. DDA simulations suggest that, like impedance-matched spheres, any impedancematched particle with a four-fold rotational symmetry also has the feature of zero backscatter, as long as the axis of symmetry is parallel to the direction of the incident radiation. This conclusion can also be obtained from theoretical analysis. Like an impedance-matched and absorbing spherical shell, an impedance-matched and absorbing aspherical shell can also be used as an invisible cloak for monostatic radar/lidar detection, provided that the above symmetry condition is satisfied. The backscatter cloak for aspherical particles, however, has less generality compared with that for spheres, as here a specific orientation of the particle is required. The DDA formalism also shows great potential to simulate the near field of a cloaked object.
This research is supported by the Office of Naval Research under contracts N00014-02-1-0478 and N00014-06-1-0069. For these grants, Dr. George Kattawar is the principal investigator. This study is also partially supported by a grant (ATM-0239605) from the National Science Foundation Physical Meteorology Program, and a grant (NAG5-11374) from the NASA Radiation Sciences Program managed previously by Dr. Donald Anderson and now by Dr. Hal Maring. For these two grants, Dr. Ping Yang is the principal investigator.
References and links
3. S. A. Cummer, B.-I. Popa, D. Schurig, D. R. Smith, and J. B. Pendry, “Full-wave simulations of electromagnetic cloaking structures,” Phys. Rev. E 74, 036621 (2006). [CrossRef]
4. D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science 314, 977–980 (2006). [CrossRef] [PubMed]
6. A. Alù and N. Engheta, “Achieving transparency with plasmonic and metamaterial coatings,” Phys. Rev. E 72, 016623 (2005). [CrossRef]
8. E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J. 186, 705–714 (1973). [CrossRef]
9. B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. 333, 848–872 (1988). [CrossRef]
10. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11, 1491–1499 (1994). [CrossRef]
11. J. D. Jackson, Classical Electrodynamics, (John Wiley & Sons, New York, 1975).
12. A. Lakhtakia, “General theory of the Purcell-Pennypacker scattering approach and its extension to bianisotropic scatterers,” Astrophys. J. 394, 494–499 (1992). [CrossRef]
13. 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]