1. Introduction

Large scale computational electromagnetic problems involving propagation through inhomogeneous media are often computationally very expensive. One of the major challenges is to reduce the computational cost of using methods like finite element or finite difference, in electrically large problems, without sacrificing the accuracy of the result. One approach for reducing the computational complexity is to downgrade the optics abstraction level from a vector EM model to a scalar optics (wave optics) model or even further down to ray tracing model. However, in some cases these approaches are unacceptable. The problem difficulty increases when complete 3D vector field results are required in real-time. Some of the motivating applications are: forward model calculation in computational imaging systems with an array of axi-symmetric lenses [1–3], radome design and distortion mitigation for radio telescopes and airborne vehicular antennas [4,5], specialized coaxial waveguide or conical feed antenna design [6,7], magnetic coil design or Fresnel plate design for wireless power transfer or super-lensing applications [8,9] and recently, in plasmonic applications such as excitonic laser [10,11] and metallic nanosphere design for plasmonic enhancement at optical frequencies [12–14]. The above mentioned electromagnetic problems are computationally hard due to the fact that the scattering structures involved are large compared to the operating wavelength (> 100λ), or a design optimization process requires repeated calculation of the scattering fields. However, in some cases the scattering structures might possess mirror/reflection, translational or rotational symmetry which can be exploited to reduce the computational cost [15]. For example, a simple circular cross-section waveguide possesses all of the above symmetries.

In this article, we focus on exploiting rotational/axial symmetry in scattering structures (which are bodies of revolution), as this symmetry reduces the dimensionality of the problem. In such structures the solutions (modes) can be represented by products of separate functions of radial, azimuthal and longitudinal variables. This in turn helps to reduce a 3D problem (ρ, ϕ, z) to a 2D problem (ρ, z) and thereby reduces the computational effort [16]. (The excitation source need not be axi-symmetric, to exploit this symmetry.) This happens due to the fact that each cylindrical harmonic propagates independently. Often, the fields of the arbitrary source are decomposed into a series of plane wave components using Fourier transform. The response from each of the obliquely incident (w.r.t z) plane wave is simulated in 2D (ρ, z) and these responses are revolved and summed to get the complete 3D solution (2.5D technique). The spectral decomposition happens in plane wave basis (Cartesian) and the evaluation of scattered fields happens in cylindrical basis.

Instead of performing spectral decomposition of the excitation source in plane wave basis, we wish to perform the decomposition in cylindrical harmonic basis (also a complete basis). The choice of cylindrical basis avoids the unnecessary coordinate transformation (which, in some cases, requires regridding or interpolation) that might be needed if plane-wave basis is used. We consider both a scalar and vector cylindrical-harmonic basis decomposition to accommodate the cases when either: the normal components of the electric and magnetic fields are known, or the transverse components of the electric field are known. We also formulate a discrete cylindrical harmonic decomposition using generalized discrete Fourier-Hankel transforms. This decomposition approach also provides a significant computational advantage in terms of constructing the overall response transfer matrix for the scattering structure involved. This sparse block diagonal transfer matrix is independent of the arbitrary source of excitation, provided, one pre-calculates the scattering response for all possible propagating modes with unit mode coefficients for a given space-bandwidth product. Our approach provides the most general methodology to exploit axial symmetry, however, some very specific numerical methods exist for cylindrical waveguides [17] and cylindrical rod scatterers [18].

In this article, we first formalize cylindrical harmonic decomposition of arbitrary fields (including both TM^z and TE^z polarizations, which are in further sections referred to as TM and TE respectively) by deriving the continuous integral equations for mode coefficients using appropriate orthogonality conditions. In Section 3, we discretize these continuous integrals using discrete Fourier-Hankel transforms (DFHT). For the discrete vector basis we formulate generalized DFHT. We also develop a composite indexing approach which facilitates the expression of these transforms in matrix form, contributing to the efficiency of the cylindrical decomposition process. In Section 4, we describe the 2.5D cylindrical harmonic simulations using the COMSOL RF module, whose scattered field results are used to compute the transfer matrix for a given scattering medium. We show four numerical validation results involving axi-symmetric scattering media namely: (a) gradient refractive index (GRIN) hemispherical Maxwell’s fish-eye lens, (b) Coordinate transformed spherical invisibility cloak with anisotropic material properties illuminated with offset vector Gaussian beam, (c) A simple ∼ 40λ diameter thin-aspheric-plano-convex lens illuminated with Gaussian vector beam, and (d) An isotropic negative index perfect lens with arbitrary excitation field. In each case, we performed a series of axi-symmetric 2D finite element simulations for all possible propagating modes to construct a transfer matrix of the scattering medium. Then this transfer matrix is multiplied by the cylindrical harmonic coefficients of the excitation field to derive the output mode coefficients at a prescribed plane. The results of this method are compared and agree well with analytical results. For the GRIN lens case (a), we also perform complete conventional 3D finite element simulations to test accuracy and compare the computational resource requirements. For the other numerical validation cases, there exist straightforward analytical results for comparison. The choice of using a cylindrical harmonic propagator in such problems is also well justified, and we compare the computational requirement of this method with an another approach, namely, the plane wave propagator. Even though we restrict our discussions to EM problems, some of the techniques and methods described in this article can also be extended to other areas of computational physics problems such as: acoustics, thermodynamics, and fluid dynamics.

2. Cylindrical harmonic decomposition

Cylindrical harmonic decomposition (CHD) is the process of decomposing an arbitrary function into basic patterns/harmonics that have simple radial and angular structures. Cylindrical harmonics, just like plane waves or spherical harmonics, form a complete orthonormal basis set. Any arbitrary field that is a solution to Helmholtz’s equation can be represented as an infinite weighted sum of these cylindrical harmonics. Often, for most of the real world signals of interest, the cylindrical harmonic expansion converges rapidly leading to good representation accuracy using a finite and relatively small number of cylindrical harmonic modes. Any arbitrary vector field can be decomposed into cylindrical harmonics as follows. Consider an arbitrary electromagnetic vector field including both TM and TE polarizations propagating in +z direction. We decompose either the z-components of both the electric and magnetic fields, or the transverse component of electric field, into cylindrical harmonic modes. We formulate the expansions so that both types use the same set of expansion coefficients. These expansion coefficients can then be propagated (from a plane z = z₀ to another plane z = z₁) by a transfer matrix and recomposed into the spatial field representation as follows.

For decomposing the field z-components we use the usual scalar basis functions employed in the analysis of cylindrical wave-guides, but the transverse electric field decomposition requires a less well-known set of vector basis functions. Below we show both the continuous and discretized versions of this change of basis. For the discrete case of the change to the scalar basis we use a Discrete Fourier-Hankel Transform (DFHT). For the discrete case of the change to the vector basis we use a generalization of the DFHT.

In what follows we assume a maximum radius over which the fields are non-negligible is given by R, and a maximum radial phase constant (or maximum radial spatial bandwidth) is given by B. The latter, might be for example, the maximum radial phase constant which permits propagating fields, i.e. the plane wave propagation constant in the medium.

2.1. Scalar basis

The expansion of arbitrary z-components of the fields at z = z₀ is given by [19],

The basis functions and the normalizing coefficients are defined separately for convenience in representing the discrete transforms below. Within each mode type (TM or TE) the basis functions are ortho-normal,

Using the orthogonality condition and the expansions above, we find the expressions for the cylindrical harmonic mode coefficients,

2.2. Vector basis

We expand the transverse electric field,

Using this condition and the field expansion above, we find the expression for the expansion coefficients as,

3. Discrete Fourier-Hankel transform

The discretized form of transforming to the cylindrical harmonic basis from the position basis requires a combined discrete Fourier transform of the angle, ϕ, and a discrete Hankel transform of the radial variable, ρ. This is seldom done because, on the surface, the Discrete Fourier-Hankel Transform (DFHT) appears to be more complex and confusing than the applying the Discrete Fourier Transform (DFT) in the Cartesian basis. A relatively simple, but useful, scheme for understanding and implementing the DFHT incorporates three composite indices: one for the cylindrical harmonic components, l, one for the input samples, i, and one for the output samples, α. Here we use bold face to denote these composite indices.

We wish to compute a set of cylindrical harmonic coefficients that span a specified range of in-plane spatial bandwidth, B. For a given domain radius, R, this implies a specific space-bandwidth product, BR. Our composite index should include all modes that lie in this space-bandwidth product, and thus satisfy β_lmn R = χ_lmn ≤ BR.

We can start from an initial set of azimuthal and radial mode indices, m and n, with maximum values given by,

The best accuracy for cylindrical harmonic coefficients results by using the roots of the radial basis functions (i.e. the first-kind Bessel functions) for selecting the radial sample points. These samples are not uniformly spaced and their values depend on both the mode type, l, and azimuthal mode number, m, of the cylindrical harmonic coefficient being computed, that is to say that the computation of each cylindrical harmonic coefficient has its own unique set of radial sample points. The density of the samples points in the spatial domain is controlled by the desired spatial bandwidth. Just as with the DFT, a certain amount of over-sampling provides better accuracy. We define the sampling bandwidth, B_s, to be larger than the spatial bandwidth of the desired cylindrical harmonics, B_s > B.

The index ranges for the azimuthal angle and radial samples are controlled by this sampling bandwidth,

To form a proper inverse, the output sample coordinates, that result from the inverse DFHT, should be the same as the input sample coordinates. However, this leads to distinct sets of sample coordinates for each cylindrical harmonic mode. Typically, one wants the total field, due to all modes, on some convenient grid of points. One could perform the inverse DFHT and then interpolate the fields from the different modes onto a single grid, or one could accomplish both operations together, by choosing a set of output sample coordinates that are mode independent, and in effect using the basis functions for the interpolation. Though any suitable grid can be used, here we show a Cartesian grid,

An example of such composite indices are shown in Table 1 on page 9, with unit domain radius, spatial bandwidth of B = 2π, sampling bandwidth of B_s = 1.1B, and grid spacing of δ = 0.5. The representative form of the transfer matrix T for this example is shown in Eq. (46) in Appendix C. We also present a pseudo-algorithm for one of the composite indicies, l, in Appendix C.

Table 1. Composite indices l, i, and α constructed for an example case: R = 1, B = 2π, B_s = 1.1B and δ = 0.5. Highlighted (grey colored rows) in l index indicate that the harmonics have space-bandwidth product less than or equal to RB. Highlighted (grey colored rows) in α index indicate that the spatial coordinate is within the domain radius R.

View Table | View all tables in this article

We can use these composite indices to form matrices for efficient computation of DFHT as;

3.1. Scalar basis DFHT

For the case of the scalar basis decomposition, the matrix elements for: the input samples, the forward DFHT, the inverse DFHT, and the output samples are given by,

3.2. Vector basis DFHT

For the case of the vector basis decomposition, the matrix elements for: the input samples, the forward DFHT, the inverse DFHT, and the output samples are given by,

One might assume, that since the l = 2 scalar basis functions are non-zero (but with zero radial derivative) at the the domain boundary, ρ = R, that one can accurately find cylindrical harmonic coefficients for H_z fields that are finite at ρ = R. In fact, the accuracy will be poor in this case. Only when fields are near zero at the boundary, can the DFHT provide accurate coefficients. The usual derivation of the discrete Hankel transform involves approximating a finite radial integration by an infinite one. This approximation is only accurate when the functions being transformed is zero outside the specified domain radius. Due to the intrinsic periodicity, an analogous problem exists for functions with non-zero amplitude at the boundary of a Fast Fourier Transform (FFT).

4. Cylindrical harmonic mode propagation (2.5D) & construction of transfer matrix T

The vector wave propagation of electromagnetic fields through an inhomogeneous medium is a well known problem and is often solved using finite difference time domain (FDTD) or finite element methods (FEM). However, this problem becomes computationally challenging when the size of this inhomogeneous structure is fairly large compared to wavelength. This forces one to exploit symmetries in the structure in order to reduce the computational load [21]. In this section we describe a technique to model axially symmetric structures using quasi-two-dimensional modeling method (2.5D modeling). This approach exploits the rotational symmetry of the structure (about the z axis) without imposing any restrictions on the source excitation field.

In section 2 and 3, we described the decomposition of arbitrary fields in terms of cylindrical harmonics. The TM and TE coefficients determined by this method can be used to completely describe the source excitation field. Similarly, the scattered fields due to the test object can again be decomposed into TM and TE cylindrical harmonic modes and their respective coefficients can be found. These scattered field coefficients form the transfer matrix T shown in Eq. (1). By combining the excitation and the scattered field coefficients appropriately with their respective basis functions one can determine the total fields due to the scattering test object. When the inhomogeneous scattering structure is rotationally symmetric (ϕ independent), a cylindrical harmonic mode of a given azimuthal mode number, m, couples only to modes of the same mode number, m. In such a structure, the electric field in cylindrical coordinate varies with the azimuthal mode m as:

 

Fig. 1 Transfer Matrix T in its sparse, block-diagonal form. The m = 0^th block is highlighted in magenta.

Download Full Size | PDF

The computational complexity of this numerical technique is mainly dependent on the total number of band-limited modes; S in a given space-bandwidth product; R × B. The total number of modes increases quadratically with space-bandwidth product and this is shown in the Fig. 2. The number of 2D cylindrical harmonic simulations and in turn the size of the transfer matrix T are directly dependent on total number of modes, S as shown in Fig. 1. An important advantage of this approach is that the 2D cylindrical harmonic simulations are independent of one another and can be completely parallelized. The computational complexity of this method can now be written as O(S).

 

Fig. 2 Shows the total number of band-limited modes; S as a function of space-bandwidth product R × B along with a quadratic fit.

Download Full Size | PDF

5. Numerical validation

The following examples are considered to numerically verify the techniques described in the previous sections.

5.1. Gradient index hemispherical Maxwell’s fish-eye lens

Maxwell’s fish-eye lens is a radially symmetric GRIN lens whose refractive index profile n, is of the form, $n(r)=\frac{2}{1+{\left(r/a\right)}^2}$, where, a is the radius of the lens and $r=\sqrt{x^2+y^2+z^2}$ is the radial coordinate. We simulate the hemispherical version of this lens excited by an offset Gaussian vector beam. For this example, we perform a detailed computational performance analysis between complete 3D full-wave and 2.5 D cylindrical harmonic propagation technique. We also compare the numerical results between the two methods. For this purpose, we chose a nominal lens radius, a = 4λ and an incident offset Gaussian full beam-waist of 2λ propagating along +z direction. (An offset Gaussian beam ensures a non-symmetric excitation field with harmonic components that have azimuthal mode numbers m greater than 1). We first simulate the blue region of the lens (one half of the 2D slice of the lens named as 2.5 D COMSOL Simulation Domain) as shown in Fig. 3 and construct the transfer matrix T based on all possible propagating modes inside this simulation domain (R=8λ). The input cylindrical harmonic mode coefficients of the source field and the transfer matrix are used to generate output plane mode coefficients. These output modes are then expanded (and also independently further propagated) using their respective basis functions to calculate the full wave vector fields post the lens. This procedure is as explained in section 4. The electric field norm on the xz plane (y = 0) is shown in Fig. 4(a). The internal fields of the lens for the 2.5D case are not shown as the transfer matrix was constructed for a plane just outside of the lens (z = 4λ) and the resulting output modes were propagated. We also simulated the full-wave complete 3D structure in COMSOL. These results are shown in Fig. 4(b). We plot the transverse plane electric field norm (xy plane) at the focal plane of the lens and at the source plane in Fig. 4(c) and (d) respectively. At the focal plane, we subtract the results from cylindrical harmonic propagator technique (2.5D) and 3D full wave simulations and show that the results are indeed accurate within a maximum error < 1 %. The normalized absolute error map at the focal plane of the lens is shown in Fig. 4(e).

 

Fig. 3 Hemispherical Maxwell’s fish-eye gradient index lens schematic for 2.5 D technique. Blue outline box indicates the simulated domain (ρ, z, ϕ = 0).

Download Full Size | PDF

 

Fig. 4 Shows the comparison between 2.5D technique and complete 3D simulations. (a). Electric field norm on xz plane (y = 0) using 2.5 D technique (Only fields pre and post the scattering medium are computed using the transfer matrix. Hence the fields inside the lens in (a) are blank). (b). Electric field norm on xz plane (y = 0) computed using conventional 3D approach in COMSOL. (c). Electric field norm at the focal plane of the lens. The white dashed line indicates the focal line. (d). Offset Gaussian beam with a beam waist=2λ propagating along +z direction is the source excitation field in both 2.5D and 3D simulations. (e). Normalized absolute error between 2.5D technique and 3D simulations at the focal plane. Electric fields are in the units of V/m.

Download Full Size | PDF

5.1.1. 3D versus 2.5D based transfer matrix

To compare a complete 3D full-wave simulation with a cylindrical harmonic propagator, using the 2.5D method and a transfer matrix approach, an example with a specific space-bandwidth product was analyzed. Performance analysis comparison between the two methods is presented in Table 2 on page 15. The cylindrical harmonic technique has superior performance in terms of computational burden, simulation time, memory usage, computational scaling, result re-usability, and parallelizability.

Table 2. Comparison of Computational Performance between complete 3D simulation and Cylindrical Harmonic Mode Propagation technique (2.5D technique) for a domain radius R = 8λ.

View Table | View all tables in this article

5.1.2. Plane wave versus cylindrical harmonic basis

As stressed in previous sections, construction of a reusable transfer matrix for a given axi-symmetric medium is one of the main goals of this article. The transfer matrix, T, for a given medium can be constructed using any complete basis, including plane waves, which comprise the most commonly used basis. Since plane waves cannot be simply represented in cylindrical coordinates, one is at first tempted to use a full-wave 3D solver (rather than an axi-symmetric solver). Such an approach is indeed more straightforward, but does not fully exploit the axial symmetry of the mediums of interest in this article. One can employ one mirror plane symmetry but the solution domains will remain three dimensional, with the associate unfavorable size scaling. Somewhat surprisingly, it is more advantageous to expand a plane wave in a cylindrical-harmonic basis and use the axi-symmetric solver as we now describe. The z-components of a plane wave incident at an arbitrary angle, θ_i (angle between the wave vector and ẑ), can be expanded in cylindrical coordinates using cylindrical harmonics as [19,21],

In order to construct the transfer matrix T for an axi-symmetric medium based on plane wave propagator, each cylindrical harmonic that contributes to a plane wave (every plane wave has a discrete value of θ_i in the range 0 to π/2) has to be propagated using an axi-symmetric solver. The 2D scattered fields due to a plane wave at z = z₁ plane is calculated by revolving the 1D fields in (ϕ = 0, z = z₁) cut-line (i.e multiplying by the factor e^jmϕ, with ϕ discretely varying from 0 to 2π) and summing all the 2D field contributions from each cylindrical harmonic. These summed 2D scattered fields due to a single plane wave input need to be regridded and interpolated to transform from (ρ, ϕ) cylindrical coordinates to (x,y) Cartesian coordinates. Then, 2D FFT can be applied on these fields to calculate the output harmonic coefficients for a particular input harmonic. (Parity conditions exist and can be exploited for −ve azimuthal mode numbers). Once the output coefficients are known, inverse FFT can be applied to calculate the vector field results post the medium. In this approach as one can very evidently note, there is a back and forth change of basis. Hence, this approach to construct T may not be the best choice in terms of computational load for a given space-bandwidth product inspite of superior, well known FFT algorithms. Table 3 illustrates the case for cylindrical harmonic propagator by comparing the computational load required to construct T using the two methods for a specific (Maxwell’s Fish-eye with R = 8λ) example.

Table 3. Comparison of computational load between Plane wave propagator and Cylindrical harmonic propagator using the same 2D axi-symmetric solver (2.5D technique) for a domain radius R = 8λ.

View Table | View all tables in this article

5.2. Spherical invisibility cloak

A spherical invisibility cloak is a transformation optics designed media whose material properties are anisotropic but yet radially symmetric. This example was chosen to show the versatility of cylindrical harmonic propagator technique in exploiting axial symmetry and in dealing with such sophisticated medium. The second-order electric permittivity tensor, ∊ (or the magnetic permeability, μ) for a spherical invisibility cloak in spherical coordinates (r, θ, ϕ) is given by [23],

 

Fig. 5 Spherical invisibility cloak schematic for 2.5 D technique. Blue outline box indicates the simulated domain.

Download Full Size | PDF

 

Fig. 6 Shows the real part of the field (in V/m) due to incident offset Gaussian vector beam on xz plane (y = 0) inside and outside the spherical cloak simulated using 2.5D technique. The subfigures below show the individual component of the material property tensors required inside the cloaking region.

Download Full Size | PDF

5.3. Free space impedance matched thin aspheric plano-convex lens

In this case, we consider a simple aspheric, thin plano-convex lens with a focal length, f = 15λ, refractive index of n = 3.58 and relative impedance Z = 1 (∊_r = μ_r = n) as the homogeneous linear scattering medium. The lens is axially symmetric with a maximum thickness of 0.8λ and a diameter of 18λ. The lens is illuminated normally by a ŷ-polarized Gaussian beam of waist diameter = 5λ at the source plane (40λ behind the lens). This example was chosen to demonstrate the versatility of this method when the transverse size and range are large compared to λ. Figure 7 shows the complete schematic diagram of the setup. The blue colored rectangular box in the figure shows the COMSOL 2.5 D simulation domain along with the boundary conditions that were applied. The simulation domain has a transverse radius R_domain = 30λ. Figure 8(a) shows the numerically computed electric field norm pre- and post-lens. The input Gaussian beam is focused at z = 15λ in front of the lens as per the design. Since a reflection-less lens in general converts a normally incident input Gaussian beam to another output Gaussian beam, one can derive the complete analytical solution for this setup using Gaussian beam optics theory [24]. Using the parameters of the considered lens and the excitation Gaussian beam, the waist radius of the output beam as a function of z was analytically calculated. Figure 8(b) shows the comparison between the numerical technique and the exact analytical result. The numerical and the analytical output beam radius are in good agreement with each other within an RMS error of 3%.

 

Fig. 7 Aspheric thin plano-convex lens schematic for 2.5 D technique. Blue outline box indicates the simulated domain.

Download Full Size | PDF

 

Fig. 8 (a). Numerically computed fields norm of electric field pre- and post-lens (in V/m). (b). Comparison between the output Gaussian beam radius as a function of z calculated numerically (black solid) and analytically (dashed red).

Download Full Size | PDF

5.4. Negative refractive index perfect lens

In this example, we choose a perfect negative index lens excited by a completely arbitrary source excitation field. A perfect lens has a straightforward analytical relationship between the image and source fields. A perfect lens (∊_r = μ_r = n = −1) acts as a focusing slab in which the image is the exact reproduction of the source (when there is no loss) [25]. The perfect lens disk has a thickness z₂ = 3.6λ and a radius R_domain = 30λ. The distances from the edge of the lens slab to the object and image plane are z₁ = 0.45λ and z₃ = 3.15λ respectively, such that z₂ = z₁ + z₃. Figure 9 shows the schematic diagram of the setup along with the COMSOL simulation domain and the boundary conditions used.

 

Fig. 9 Perfect Lens schematic for 2.5D technique. Blue outline box shows the simulated domain.

Download Full Size | PDF

Figure 10(a) shows the source plane electric field norm of the excitation field. Since the image is the object for a perfect lens with no loss, the fields at the image plane should be exactly equal to the source fields. Figure 10(b) shows the difference between the image and source fields and the absolute peak error about 1.5%.

 

Fig. 10 (a) Shows the source plane electric field norm (in V/m) of the chosen arbitrary excitation source at z = 0 plane. (b) Shows the difference map between the image and the object plane electric field norms (in V/m).

Download Full Size | PDF

6. Conclusion

In this article we develop a complete formulation for both continuous and discrete cylindrical harmonic decomposition. We formulate this decomposition method for both scalar and vector basis fields. Such harmonic decomposition methods (like Fourier transforms and spherical harmonics) play a vital role in many computational physics problems. Some future developments regarding the decomposition will include reducing the computational burden further and parallelizing the decomposition technique on a Graphical Processing Unit (GPU). We use this decomposition method to propagate the individual harmonics through an axi-symmetric medium using a 2D axi-symmetric EM solver as a part of a pre-computational step (2.5D technique). The fields from this solver are inturn used to construct an efficient, modular, and re-useable transfer matrix, which can be later used to calculate the total scattered fields due to the inhomogeneous axi-symmetric medium for any given arbitrary incident excitation and a finite space-bandwidth product. We also present a composite index ordering strategy to perform the required forward and/or inverse DFHT. This ordering helps to reduce the operations to simple matrix multiplications. This overall approach was numerically verified on well-known and relevant axi-symmetric media. We also compare our method with other possible techniques in terms of computational load. The exploitation of axial symmetry along with transfer matrix approach using CHD help to significantly reduce the computational burden/resources without sacrificing accuracy.

A. Orthogonality Properties of Bessel and Exponential Functions

The following orthogonality properties of the Bessel and exponential functions are used in the article,

(43)$$\underset{0}{\overset{R}{\int }}{J}_m\left({\chi }_{1mn}\hspace{0.17em}\hspace{0.17em}\rho \right)\hspace{0.17em}\hspace{0.17em}{J}_m({\chi }_{1m{n}^{\prime }}\hspace{0.17em}\rho )\rho d\rho =\frac{R^2}{2}{J}_{m+1}^2({\chi }_{1mn}){\delta }_{n{n}^{\prime }}$$

(44)$$\underset{0}{\overset{R}{\int }}{J}_m\left({\chi }_{2mn}^{\prime }\rho \right)J_m\left({\chi }_{2m{n}^{\prime }}^{\prime }\rho \right)\rho d\rho =\frac{R^2}{2}\left[J_m^2\left({\chi }_{2mn}^{\prime }\right)-J_{m-1}\left({\chi }_{2mn}^{\prime }\right)J_{m+1}\left({\chi }_{2mn}^{\prime }\right)\right]{\delta }_{n{n}^{\prime }}$$

(45)$$\underset{-\pi }{\overset{\pi }{\int }}\text{exp}(-jm\varphi )\text{exp}(j{m}^{\prime }\varphi )d\varphi =2\pi {\delta }_{m{m}^{\prime }}$$

B. Parity Condition between positive and negative azimuthal mode numbers

The following Tables 4 and 5 show the parity condition between positive and negative azimuthal mode numbered basis functions for TM and TE modes.

Table 4. Parity Condition for TM modes.

View Table | View all tables in this article

Table 5. Parity Condition for TE modes.

View Table | View all tables in this article

C. Composite Indices for DFHT

Based on the example composite indices shown in Table 1 on page 9, the transfer matrix T takes the representative form shown in Eq. 46, where, ⋈ and # in T matrix indicates the TM (l′ = 1) and TE (l′ = 2) mode coefficients respectively. The rest of the values in the 17 × 17 transfer matrix are zeros, making it a sparse, block-diagonal matrix.

The functional form of the composite indices [l, i, α] mapping are neither illuminating nor easy to specify in mathematical notation. It is best to generate a vector of triplets that store the mapping and use elements of this vector for building up the needed matrices. Some pseudo-code for generating this composite index vector l is shown in Algorithm 1.

(46)$$\left[\mathbf{b}\right]={\left[\begin{array}{ccccccccccccccccc}\#& & & & & & & & & & & & & & & & \\ & \#& & & & & & & & & & & & & & & \\ & & \bowtie & \bowtie & & & & & & & & & & & & & \\ & & \#& \#& & & & & & & & & & & & & \\ & & & & \bowtie & \bowtie & \bowtie & & & & & & & & & & \\ & & & & \#& \#& \#& & & & & & & & & & \\ & & & & \#& \#& \#& & & & & & & & & & \\ & & & & & & & \bowtie & \bowtie & \bowtie & \bowtie & & & & & & \\ & & & & & & & \#& \#& \#& \#& & & & & & \\ & & & & & & & \#& \#& \#& \#& & & & & & \\ & & & & & & & \#& \#& \#& \#& & & & & & \\ & & & & & & & & & & & \bowtie & \bowtie & \bowtie & & & \\ & & & & & & & & & & & \#& \#& \#& & & \\ & & & & & & & & & & & \#& \#& \#& & & \\ & & & & & & & & & & & & & & \bowtie & \bowtie & \\ & & & & & & & & & & & & & & \#& \#& \\ & & & & & & & & & & & & & & & & \#\end{array}\right]}_{17\times 17}{\left[\begin{array}{c}{a}_2,-4,1\\ a_2,-3,1\\ a_1,-2,1\\ a_2,-2,1\\ a_1,-1,1\\ a_2,-1,1\\ a_2,-1,2\\ a_{1,0,1}\\ a_{1,0,2}\\ a_{2,0,1}\\ a_{2,0,2}\\ a_{1,1,1}\\ a_{2,1,1}\\ a_{2,1,2}\\ a_{1,2,1}\\ a_{2,2,1}\\ a_{2,3,1}\end{array}\right]}_{17\times 1}$$

Algorithm 1. Finding Composite Index l

View Table | View all tables in this article

Funding

This work was supported by the Department of Homeland Security, Science and Technology Directorate (Contract No. HSHQDC-12-C-00049). The published material represents the position of the authors and not necessarily that of the DHS or S&T.

Acknowledgments

The support and resources from the Center for High Performance Computing at the University of Utah are gratefully acknowledged. Authors would like to acknowledge Richard Edwards for some useful comments on the manuscript.

References and links

1. O. Furxhi, D. L. Marks, and D. J. Brady, “Echelle crossed grating millimeter wave beam scanner,” Opt. Express 22, 16393–16407 (2014). [CrossRef]   [PubMed]  

2. O. S. Cossairt, D. Miau, and S. K. Nayar, “Scaling law for computational imaging using spherical optics,” J. Opt. Soc. Am. A 28, 2540–2553 (2011). [CrossRef]  

3. S. Venkatesh, N. Viswanathan, and D. Schurig, “W-band sparse synthetic aperture for computational imaging,” Opt. Express 24, 8317–8331 (2016). [CrossRef]   [PubMed]  

4. M. Meeks and J. Ruze, “Evaluation of the haystack antenna and radome,” Antennas and Propagation, IEEE Transactions on 19, 723–728 (1971). [CrossRef]  

5. G. Crone, A. Rudge, and G. Taylor, “Design and performance of airborne radomes: A review,” Communications, Radar and Signal Processing, IEE Proceedings F 128, 451–464 (1981). [CrossRef]  

6. M. Ibanescu, Y. Fink, S. Fan, E. L. Thomas, and J. D. Joannopoulos, “An All-Dielectric coaxial waveguide,” Science 289, 415–419 (2000). [CrossRef]   [PubMed]  

7. B. Thomas, G. L. James, and K. J. Greene, “Design of wide-band corrugated conical horns for cassegrain antennas,” Antennas and Propagation, IEEE Transactions on 34, 750–757 (1986). [CrossRef]  

8. G. Lipworth, J. Ensworth, K. Seetharam, D. Huang, J. S. Lee, P. Schmalenberg, T. Nomura, M. S. Reynolds, D. R. Smith, and Y. Urzhumov, “Magnetic Metamaterial Superlens for Increased Range Wireless Power Transfer,” Scientific Reports 4, 3642 (2014). [CrossRef]   [PubMed]  

9. R. Merlin, “Radiationless electromagnetic interference: Evanescent-Field lenses and perfect focusing,” Science 317, 927–929 (2007). [CrossRef]   [PubMed]  

10. Y. Ye, Z. J. Wong, X. Lu, X. Ni, H. Zhu, X. Chen, Y. Wang, and X. Zhang, “Monolayer excitonic laser,” Nat Photon 9, 733 (2015). [CrossRef]  

11. Y. Zhao, J. Wyrick, F. D. Natterer, J. F. Rodriguez-Nieva, C. Lewandowski, K. Watanabe, T. Taniguchi, L. S. Levitov, N. B. Zhitenev, and J. A. Stroscio, “Creating and probing electron whispering-gallery modes in graphene,” Science 348, 672–675 (2015). [CrossRef]   [PubMed]  

12. J. B. Pendry, A. Aubry, D. R. Smith, and S. A. Maier, “Transformation optics and subwavelength control of light,” Science 337, 549–552 (2012). [CrossRef]   [PubMed]  

13. C. Ciraci, R. T. Hill, J. J. Mock, Y. Urzhumov, A. I. Fernández-Dominguez, S. A. Maier, J. B. Pendry, A. Chilkoti, and D. R. Smith, “Probing the Ultimate Limits of Plasmonic Enhancement,” Science 337, 1072–1074 (2012). [CrossRef]   [PubMed]  

14. Y. Ou, D. Pardo, and Y. Chen, “Fourier finite element modeling of light emission in waveguides: 2.5-dimensional fem approach,” Opt. Express 23, 30259–30269 (2015). [CrossRef]   [PubMed]  

15. C. C. Johnson, Field and Wave Electrodynamics (McGraw-Hill, 1965).

16. A. F. Peterson, S. L. Ray, and R. Mittra, Computational Methods for Electromagnetics, Volume 4 of IEEE Press Series on Electromagnetic Wave Theory (Wiley, 1998).

17. R. Song, J. Zhu, and X. Zhang, “Full-vectorial modal analysis for circular optical waveguides based on the multidomain chebyshev pseudospectral method,” J. Opt. Soc. Am. B 27, 1722–1730 (2010). [CrossRef]  

18. S. Feng, K. Halterman, P. L. Overfelt, and D. Bowling, “Cyclic sommerfeld resonances in nanorods at grazing incidences,” Opt. Express 17, 19823–19841 (2009). [CrossRef]   [PubMed]  

19. C. A. Balanis, Advanced Engineering Electromagnetics, 2nd Edition (Wiley, 2012).

20. S. Ali, W. Chew, and J. Kong, “Vector hankel transform analysis of annular-ring microstrip antenna,” IEEE Transactions on Antennas and Propagation 30, 637–644 (1982). [CrossRef]  

21. C. Ciraci, Y. Urzhumov, and D. R. Smith, “Far-field analysis of axially symmetric three-dimensional directional cloaks,” Opt. Express 21, 9397–9406 (2013). [CrossRef]   [PubMed]  

22. COMSOL Multiphysics, RF Module User’s Guide (COMSOL, 2012).

23. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science 312, 1780–1782 (2006). [CrossRef]   [PubMed]  

24. B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics (Wiley, 2007).

25. J. B. Pendry, “Negative Refraction Makes a Perfect Lens,” Physical Review Letters 85, 3966–3969 (2000). [CrossRef]   [PubMed]