1. Introduction

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

where

are the tangential components, and

are the radial components of the permittivity and permeability, respectively, with ε₀ and µ₀ the free space permittivity and permeability, respectively, and R₁ and R₂ the inner and outer radii of the spherical cloak, respectively. Ray-tracing simulations [1, 2] supporting the conclusions by Pendry et al. [1] were reported in the geometric optics limit. Full-wave finite-element simulations [3] 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 [4] using artificially structured metamaterials at microwave frequencies. A rigorous solution to Maxwell’s equations for a cloaked sphere has been reported [5] 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. [6] 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. [7] noticed that the cloaking condition is not very strict if only backscatter is concerned. They showed that when the impedance $Z=\sqrt{\frac{\mu }{\epsilon }}$ was equal to the impedance of the surrounding medium, Z = Z₀, 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.[7] 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 [11] 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 [12], 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 P_j=α̿_j E _ext,j, where P_j 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 H_{ext, j}, where M_j is the instantaneous magnetic dipole moment and H_{ext, 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 [11] 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) [13] 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 E₀ × H₀ = k, the wave vector of the incident radiation. In the second term, $-{\mathbf{A}}_{jk}^{(ee)}{\mathbf{P}}_k$ and $-{\mathbf{A}}_{jk}^{(eh)}{\mathbf{M}}_k$ 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 $-{\mathbf{A}}_{jk}^{(he)}{\mathbf{P}}_k$ and $-{\mathbf{A}}_{jk}^{(hh)}{\mathbf{M}}_k$ are the contributions to the magnetic field at position j due to the electric and magnetic dipoles at position k, given by [11] (pp. 413)

where r_jk = r_j − r _k and r_jk = |r_jk|. Notice here that the electric dipole moments and magnetic dipole moments are coupled with each other.

Let A_jk 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, A_jk, is a 3×3 matrix and that $-{\mathbf{A}}_{jk}^{(ee)}$ and $-{\mathbf{A}}_{jk}^{(hh)}$ are symmetric, while $-{\mathbf{A}}_{jk}^{(eh)}$ and $-{\mathbf{A}}_{jk}^{(he)}$ 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 M_j’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.

 

Fig. 1. Operations on a four-fold symmetric particle. 1, 2, 3, and 4 are a group of points that rotate into each other. E and H are the incident electric and magnetic fields. The operation from (I) to (II) is a π/2 rotation of the coordinate system; the operation from (II) to (III) is a relabeling of E and H; (III) and (I) are identical except for a rotation of the four points. The coordinate system is as shown at the lower left corner, with the z-axis pointing out from the paper.

Download Full Size | PDF

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

P_1′ = (−y,x,z), M_1′ = (−ν,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

$\sum _{k=1}^N\left[\begin{array}{cc}{\mathbf{A}}_{\mathrm{jk}}^{\left(\mathrm{ee}\right)}& {\mathbf{A}}_{\mathrm{jk}}^{\left(\mathrm{eh}\right)}\\ {\mathbf{A}}_{\mathrm{jk}}^{\left(\mathrm{he}\right)}& {\mathbf{A}}_{\mathrm{jk}}^{\left(\mathrm{hh}\right)}\end{array}\right]·\left[\begin{array}{c}\mathbf{P}{\prime }_k\\ \mathbf{M}{\prime }_k\end{array}\right]=\left[\begin{array}{c}\mathbf{H}{″}_{\mathrm{inc},j}\\ -\mathbf{E}{″}_{\mathrm{inc},j}\end{array}\right],$

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

Now compare Fig. 1(III) with Fig. 1(I). The only difference is that point 2 is now labeled as 1″. By this means, we can determine the local electric and magnetic moments at 2 in Fig. 1(I), which are

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 [11]

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 E₁₂₃₄(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 E₁₂₃₄’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 [7] for coated spheres with matched impedances. Here for both layers the refractive index is $n=\sqrt{\epsilon \mu }$ , and the impedance is $Z=\sqrt{\frac{\mu }{\epsilon }}$ . In our comparison, we used the following parameters: the inner core is a regular dielectric medium with n₁ = 1.5 and Z₁ = 1/1.5, or $\sqrt{\epsilon }=1.5$ and $\sqrt{\mu }=1$ ; the outer coating is an absorbing medium with unit impedance, n₂ = 2+i and Z₂ = 1, or $\sqrt{\epsilon }=\sqrt{\mu }=\sqrt{2+i}$ , with $i=\sqrt{-1}$ being the imaginary unit. The size parameters for the inner and outer layers are denoted by x=2πR ₁/λ and y=2πR ₂/λ, respectively. An inter-dipole spacing d~λ/40 was used in the DDA simulations. The resulting P₁₁, 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, P₁₁ 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.

 

Fig. 2. Comparison of P₁₁ for coated spheres computed from DDA and Mie methods. Two cases, x = 2.5, y = 5, and x = 4, y = 8 were simulated. For the inner core, n₁ = 1.5 and Z₁ = 1/1.5, for the outer coat, n₂ = 2+i and Z₂ = 1.

Download Full Size | PDF

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 a_eff = (3V/4π)^1/3 with V the volume of the aspherical particle. Shown in Fig. 3 is P₁₁ 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 $\sqrt{\epsilon }=\sqrt{\mu }=\sqrt{1.5}$ . For comparison, results for a regular dielectric particle (n = 1.5 and Z = 1/1.5, or $\sqrt{\epsilon }=1.5$ and $\sqrt{\mu }=1$ ) 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 P₁₁ in the two scattering planes are appreciably different. It turned out that for an impedance-matched cubic particle in this specific orientation, the simulated P₁₁ 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.

 

Fig. 3. P ₁₁ for cubic particles with size parameter x = 8 and n = 1.5, Z = 1. Results for a particle with the same shape but n = 1.5, Z = 1/1.5 are also included for comparison.

Download Full Size | PDF

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 P₁₁’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 (n₁ = 1.1 and 1.5, Z₁ = 1/n₁), and the absorbing outer cubic coating has a size parameter y=8 and a unit impedance (n₂ = 2+i and Z₂ = 1). The results are shown in Fig. 5, where the blue curve corresponds to the n₁ = 1.1 case, and the red curve corresponds to the n₁ = 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⁻⁷, 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 [1]. 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 ₁/R₂ is close to 1. For example, R₁ = 0.3R₂ requires n_t = 1.43, which is a moderate value; while R₁ = 0.5R₂ requires n_t = 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.

 

Fig. 4. (a) The projections on the x-y plane of various homogeneous impedance-matched columns, an octagonal column, a hexagonal column, a cube, and a rectangular column, oriented with one base facing the incident radiation. (b) P ₁₁ in the ϕ = 0° plane for the above columns. The refractive index was n = 1.21. Inset shows the same results from 160° to 180°.

Download Full Size | PDF

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 R₂ = 8λ/2π, and the results for two cases, R₁ = 0.3R₂ and 0.5R₂ 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 R₁ = 0.3R₂ 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 R₁=0.5R₂ 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. [3] 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 $Q_{\text{sca}}=C_{\text{sca}}/(\pi R_2^2)$ for the cases shown in Fig. 6 are on or below the order of 10^-4 (6.13×10^-5 in case (a) with R₁ = 0.3R₂ and 4.43×10^-4 in case (b) with R₁ = 0.5R₂), which are sufficiently small to exhibit the nature of invisibility. For comparison, the scattering efficiency of the uncloaked sphere in case (c) is Q_sca = C_sca/(πR²) = 3.46. We use the outer radius R₂ in the calculation of Q_sca as we are comparing cases with the same R₂. It can also be noticed that the DDA calculations gave a better simulation of the zero scattering for the cloak with R₁ = 0.3R₂, in which case n_t is smaller.

 

Fig. 5. Scattering off coated-cubes with a regular dielectric core (Z₁ = 1/n₁) and an absorbing shell with unit impedance (n₂ = 2+i and Z₂ = 1). The size parameters are x = 4 and y = 8. Results for n₁ = 1.1 and 1.5 cases are shown. The dash curves show the scattering off the dielectric cores if naked for comparison.

Download Full Size | PDF

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 P₁₁ for a cloaked sphere with R₂ = 5λ/2π, R₁ = 0.2R₂, and various refractive indices in the cloaked region. The P₁₁ for uncloaked cores are also shown for comparison. The corresponding scattering efficiencies $Q_{\text{sca}}=C_{\text{sca}}/(\pi R_1^2)$ 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 P₁₁ and Q_sca 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.

Table 1. Simulated scattering efficiencies of cloaked spheres and the corresponding uncloaked cores with various refractive indices n.

View Table

 

Fig. 6. Electric field distribution in the vicinity of cloaked and uncloaked spheres, with scales in unit of λ. (a) A cloaked sphere with R₂ = 8λ/2π and R₁ = 0.3R₂; (b) A cloaked sphere with R₂ = 8λ/2π and R₁ = 0.5R₂. (c) An uncloaked homogeneous impedancematched sphere with ε = µ = 1.21 and R = 8λ/2π. Black circles indicate the inner and outer boundaries of the cloak, and the incident radiation propagates from the left to the right.

Download Full Size | PDF

5. Conclusions

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.

 

Fig. 7. P₁₁ for a cloaked sphere with R₂ = 5λ/2π and R₁ = 0.2R₂. Refractive indices for the cloaked region n = 1.3, 1.5 and 1.7 were simulated. Dash lines are P₁₁’s for corresponding uncloaked spheres.

Download Full Size | PDF