Far-field Analysis of Axially Symmetric Three-dimensional Directional Cloaks References and Links

Axisymmetric radiating and scattering structures whose rotational invariance is broken by non-axisymmetric excitations present an important class of problems in electromagnetics. For such problems, a cylindrical wave decomposition formalism can be used to efficiently obtain numerical solutions to the full-wave frequency-domain problem. Often, the far-field, or Fraunhofer region is of particular interest in scattering cross-section and radiation pattern calculations; yet, it is usually impractical to compute full-wave solutions for this region. Here, we propose a generalization of the Stratton-Chu far-field integral adapted for 2.5D formalism. The integration over a closed, axially symmetric surface is analytically reduced to a line integral on a meridional plane. We benchmark this computational technique by comparing it with analytical Mie solutions for a plasmonic nanoparticle, and apply it to the design of a three-dimensional polarization-insensitive cloak. Modified field enhancement and extinction by plasmonic nanowire dimers due to nonlocal response, " Opt. Finite-element analysis of complex axisymmetric radiating structures, " IEEE Trans. Isotropic-medium three-dimensional cloaks for acoustic and electromagnetic waves, " J. Scattering of a plane wave from a circular dielectric cylinder at oblique incidence, " Can. A rigorous analysis of plane-transformed invisi-bility cloaks, " IEEE Trans. Cross-section comparisons of cloaks designed by transformation optical and optical conformal mapping approaches, " J. Low-loss directional cloaks without superluminal velocity or magnetic, " Opt.


Introduction
The computational problem of electromagnetic wave propagation over large distances (e.g., hundreds of wavelengths) through media with arbitrary inhomogeneity can be challenging.In the context of electromagnetic cloaks and other transformation optical media, for example, structures with specific spatial distributions of the electric permittivity and magnetic permeability tensors must often be simulated repeatedly during optimization cycles.In the context of nanophotonics, one is often interested in solving a multiscale physics problem that takes into account, for example, the microscopic features of an electromagnetic system via either classical nonlocal [1] or quantum [2] effects; in either case, the requirement of supporting length scales from the microscopic (sub-nanometer) to the electromagnetic (micron) scale places a considerable burden on the computational requirements.
When the scattering geometry to be simulated possesses underlying symmetries-such as rotational, translational or reflection-one may choose an appropriate basis in which it is possible to separate the partial differential equations governing electromagnetic waves and other contributing physics into sets of differential equations that depend on a reduced number of variables.The particularly relevant example here is that of axially symmetric geometries, where a natural basis is the set of cylindrical harmonics.In linear electrodynamics, each of the cylindrical harmonics represents an independent solution to a reduced wave equation, characterized by distinct field patterns in the ρz-plane.
For scattering problems involving axisymmetric geometries in which the excitation field breaks the axial symmetry, such as a plane wave scattering from a sphere or cylinder, a solution can be computed by expanding the excitation field in cylindrical harmonics and computing the response for each of the harmonics separately.The advantage in doing so is that instead of solving an N ρ × N z × N φ sized problem, where the indices indicate the number of mesh points in each of the dimensions, it is instead only necessary to solve m max two-dimensional problems of size N ρ ×N z .The reduction in dimensionality typically results in extremely large computational savings in terms of space and processing time.
The simulation of axisymmetric geometries by summing cylindrical harmonics has been successfully applied in numerous radiation problems.An axisymmetric radiating structure with a non-axisymmetric source represents an important class of antennas.Often, these radiating structures are electrically large, necessitating efficient and accurate computational techniques [3].Another relevant problem in electromagnetism that has been studied for several years is the modification of the radiation pattern of an antenna covered by an axisymmetric radome [4].
More recently, structures emerging in the fields of transformation optics and nanophotonics have pushed the need to exploit axisymmetric structures even further.For example, Urzhumov et al. designed a directional, three-dimensional isotropic-medium cloak by revolving a twodimensional conformal cloak [5].Likewise, Ciracì and co-workers successfully made use of the axisymmetric technique to study systems of nanoparticles strongly coupled to metallic films, in which the microscopic nature of electrons inside the metal particles was taken into account [6].
In many of these problems, the far-field, radiating fields are of primary interest.In full-wave finite-element or finite-difference simulation approaches, the sources or scattering elements must be simulated.Thus, even for very large simulation domains, it is typically impractical to achieve spatial scales beyond the near-field or radiating near-field of the source.To determine the far-field, then, typically an integral is performed over a three-dimensional surface surrounding the scattering object that propagate the near-to the far-fields.This integral convolution, however, is inefficient in that it requires the evaluation of a sum over all the harmonics for each integration point.It is possible to reduce remarkably the computational load by analytically taking out the azimuthal dependence of the solution in the far-field expression.The surface integral thus can be reduced to a line integral on the cross-section domain boundary.
In this paper, we introduce a computational approach for the evaluation of the scattered farfields corresponding to axisymmetric structures, deriving an expression analogous to the Chu-Stratton formula by expanding the far-field in terms of cylindrical harmonics.We illustrate the technique in the design of a directional invisibility cloak.

Quasi-two-dimensional modeling
We first present the main aspects of the modeling technique employed in this paper.We introduce a general formulation of quasi-two-dimensional modeling method for vector electromagnetic fields, which we refer to as 2.5D modeling.In the 2.5D formulation, all fields are decomposed in terms of the azimuthal mode number m ∈ Z.The decomposition proceeds by expanding the fields in cylindrical harmonics with respect to the variable φ : If the geometry is φ -independent (i.e. it is rotationally symmetric), each cylindrical harmonic propagates independently.This can be easily demonstrated by direct substitution of Eqs.(1) into Maxwell's equations and considering that D (m) = ε(ρ, z)E (m) and B (m) = µ(ρ, z)H (m) .For each azimuthal number m, an mth wave equation can be written as: where the operator ∇ (m) × is obtained from the curl operator by substituting derivatives with respect to φ with a factor im in cylindrical coordinates.The set of Eqs. ( 2) can be solved on the two-dimensional cross-section of the simulation domain using a standard scattered-field formulation for any arbitrary excitation.For many wave forms of interest, the cylindrical harmonic expansion converges rapidly, and therefore it can be truncated at a relatively small m = m max .An incident plane wave typically serves as a convenient field excitation.A plane wave can be decomposed in cylindrical harmonics as follows.Consider a wave polarized such that the magnetic field is transverse to the z-axis (TM z polarization), i.e.H z = 0 as shown in Fig. 1(a) and the electric field is: In the previous expression θ i is the angle between the wave vector and the z-axis.Taking advantage of the periodicity with respect of φ , it is possible to expand the z-component of E as Fourier series [7]: A similar expression could be derived for E x .However this would introduce differing factors of sin φ and cos φ in both components E φ and E ρ , violating our requirement that the azimuthal variations must be the same for all fields.It is more convenient to derive all remaining cylindrical-coordinate components directly from E z .Using Maxwell's curl equations we obtain [7]: A parity condition relating positive and negative azimuthal number exists, which further reduces the computational load by a factor of two.In general, the parity condition depends on the particular field component considered.Total fields in the three-dimensional domain can be obtained by revolution around z-axis by reintroducing the φ -dependent factor e imφ in the sums of Eqs.(1).We solve Eqs.(1) numerically using a commercially available platform based on finiteelement method (the RF module of COMSOL Multiphysics).The scattering-field formulation is used to specify the incident fields given by Eqs. ( 4), ( 5) and (6).

Far-field calculation
In this paragraph we present a general far-field expression for arbitrary axially symmetric structures.The far-field E f ar along a given direction r can be expressed as an integral of the nearfields E and H over an arbitrary surface S enclosing the scattering object using the well-known Stratton-Chu formula [8]: In the above formula n is the unit normal to S, k is the wave number and η = µ/ε is the impedance.As we show in the Appendix A, it is possible to analytically perform the integration over the azimuthal variable φ ∈ [0, 2π].The result of this integration is: where ds is the line element on the surface S corresponding to φ = constant.The linear kernels K (m) e and K (m) h are defined as: All Bessel functions are evaluated in ξ = kρ sinθ , that is J m = J m (kρ sinθ ).Details of the above derivation are included in the Appendix A. What is important here is that by explicitly assuming the azimuthal dependence, it is possible to analytically integrate along the revolved direction φ .Thus, the surface integration is reduced to a much faster line integration along the boundary of the revolved cross-section domain.This expression is valid for any arbitrary structures that possess axial symmetry.

Validation
In order to validate our derivation we consider the simple case of electromagnetic scattering from a metallic nanoparticle.Analytical solutions for spherical particles can be easily found using the Mie solutions of Maxwell's equations [9].Consider a plane wave propagating along the z-direction toward negative values, impinging on a dielectric sphere and polarized such that the electric field lies along the x-direction (equivalently, on the plane φ = 0), as shown in Fig. 1(a).From the optical theorem, the extinction cross-section σ ext is found from the far-field amplitude in the forward direction E f ar θ (0, π) [9]: The absorption cross-section is obtained by normalizing the absorbed power by the incident wave intensity, σ abs = W abs /I 0 ,, where W abs can be easily calculated from the integral: performed over the volume Ω occupied by the particle.Writing the fields using similar expressions as in (1), gives: where we denoted the revolved cross-sectional area as Ω .We note that non-zero contributions exists only for n = m, for which the integration over the azimuthal variable gives a factor of 2π.
The previous integral finally reduces to: The volume integration is thus reduced to an integration across the revolved cross-section area only.The total scattering cross-section is then obtained as σ sca = σ ext − σ abs .Figure 1(b) shows the comparison of absorption, extinction, and scattering efficiencies calculated using the standard Mie theory and those obtained by 2.5D simulations using Eq. ( 8) for the gold nanosphere of radius R = 30 nm.Notice that the two curves coincide exactly.Figure 1(c) shows the scattering differential cross-section in the plane of polarization at the electric field (E-plane) and the magnetic field (H-plane) respectively, for the wavelength corresponding to the localized surface plasmon resonance.The scattering pattern reproduces the characteristic donut shape of the dipolar radiation.

Unidirectional 3D cloak
Electromagnetic invisibility of objects large in comparison with the wavelength has been attracting enormous interest since the proposal of Pendry et al. [10].Three-dimensional cloaking devices for objects of large diameter (relative to a reference wavelength) are difficult to design in part due to the computational demands associated with solving the full-wave problem.The key properties and ingredients of both omnidirectional [10] and unidirectional [11][12][13] cloaks can be understood using geometries with axial symmetry [5].Reducing computational time per solution of the direct problem is particularly important for optimization (design) problems.The 2.5D formalism is therefore a natural instrument for the design and performance analysis of eikonal-limit cloaks [14].Here, we extend the transformation optics-based unidirectional cloak proposed in Ref. [15] to three dimensions, and provide an accurate numerical demonstration of a polarization-insensitive invisibility device.Consider the geometry depicted in Fig. 2(a).The cloak consists of three triangular regions revolved around the vertical axis z and filled by homogeneous and anisotropic materials.The double-cone shaped cloaked object of height 2a and width 2b is coated by a perfectly conducting infinitely thin film.The surrounding cloaking device has a cylindrical boundary of height 2a and width 2B = 2a.We refer to these sectors as regions I, I , and II, as shown in Fig. 2(a).Using standard transformation optics techniques [10], and applying the following piecewise-linear, continuous transformation: where χ = B/(B − b) and γ = b/a, we obtain that the non-zero material properties are given by: I : Region I is filled with a mirror imaged material from region I, and µ = ε.A full-wave simulation performed using the 2.5D technique shows that the near-field scattering from the cloaked object appears strongly reduced (see Fig. phase accumulation errors are present.Performing the far-field integration described above, we are able to determine exactly the expected reduction in scattering cross-section in the far-field.
If we define the geometrical cross-sectional area as A = πb 2 , the total scattering cross-sections are σ sca = 2.080A and σ sca = 0.077A for the uncloaked and cloaked geometries respectively.The cloaking device reduces the total cross-section to less than 4% of the original scatterer.In Fig. 3 we show the differential scattering cross-section of the uncloaked object compared to the case where the object is surrounded by the cloaking device.The cloaking reduces the differential scattering at least of one order of magnitude in all directions with respect to uncloaked object.

Conclusion
We have developed a semi-analytical formula for calculating far-field properties of arbitrary axially symmetric structures, implementing it as an extension of state-of-the-art commercial software in electromagnetic wave propagation.We tested our technique against analytical results, showing excellent accuracy.The study of a simple spherical scatterer provides ample support as a benchmark, suggesting the applicability of the method for more complex geometries, such as transformation optical structures.
We have used the 2.5D method as a tool to study the properties of a novel directional polarization-independent cloaking geometry.An analysis of scattering properties of this cloak in the far-field region shows that the cloak reduces the total scattering to an amount less than 4% that of the uncloaked object.The generality of our approach is suitable for any numerical techniques and any kind of driving fields, including off axis incidence of linearly polarized waves.

A. Appendix: Far-field formula
We present here the detailed derivation of Eq. 8 in the main text.We begin by rewriting Eq. 7 as: such that the integral is reduced to the form: where V = n × E or V = n × H.In the 2.5D formulation all the fields are expressed in the orthonormal basis C = ρ, φ , ẑ and can be decomposed in terms of the azimuthal mode number m as in Eq. ( 1).Substituting an analogous expression for V gives: where M C is the transformation matrix from C to the standard basis: −1 ] T = M C , so that the the product to the matrix M C is associative with respect to the cross product.
Performing the products r • r and C gives respectively: where we have written the vectors in a way that makes the dependence on the φ coordinate explicit.Before proceeding in our calculation it is useful to notice that since all fields have the same azimuthal dependence, e imφ , we can safely assume that also the far-field will show the same behavior.We can write then: The calculation of the far-field can be reduced to determining the far-field coefficients E (m) f ar in the half-plane corresponding to φ = 0.With this in mind and substituting the Eqs.( 21) and (22) into Eq.( 19) we find: 3 (ρ , z )e i(m+1)φ e ikz cos θ e iξ cos φ dS , where we have assumed r ∈ {φ = 0}, ξ = kρ sin θ and: In order to solve the integrals of Eq. ( 24), consider the following integral: 2π 0 e iξ cos ϕ e imϕ dϕ = 2πi m J m (ξ ), where J m is the m-th order Bessel function of the first kind.Using Eq. ( 26) we can easily perform the integrations of Eq. ( 24) over φ ∈ [0, 2π].Considering that the surface element can be written as dS = ρ dφ ds gives: where all Bessel's functions are evaluated in ξ = kρ sin θ , that is J m = J m (kρ sin θ ).The remaining integration is performed along the boundary s of the domain cross-section.The quantity in the brackets becomes then:  (33) In writing Eq. ( 31) we omitted the component E (m) r, f ar since it is equally zero.The total far-field E f ar (φ , θ ) can be found using the expression (23).This expression is valid for any arbitrary structures that possess axial symmetry.

Fig. 1 .
Fig.1.Scattering properties of a gold nanosphere of radius R = 30 nm are calculated in the far-field implementing Eq. (8).(b) Comparison between analytical calculations based on Mie solutions (solid lines) and COMSOL 2.5D simulations (markers).(c) Differential scattering cross-section for a gold nanosphere on the E-plane (blue) and H-plane (green) as defined in (a) for θ i = 0.The wavelength of the incident radiation is λ 0 = 520 nm.

(Fig. 3 .
Fig. 3. Scattering properties of the cloaking device of Fig. 2. Normalized differential scattering cross-section of the uncloaked (left) and cloaked (right) object.The top and bottom plots refer to the scattering measured at the H-plane and E-plane respectively.