Additional modes in a waveguide system of zero-index-metamaterials with defects

Zero-index-metamaterials (ZIM) have drawn much attention due to their intriguing properties and novel applications. Particularly, in a parallel plated ZIM waveguide system with defects, total reflection or transmission of wave can be achieved by adjusting the properties of defects. This effect has been explored extensively in different types of ZIM (e.g., epsilon-near-zero metamaterials, matched impedance ZIM, or anisotropic ZIM). Almost all previous literatures showed that only monopole modes are excited inside the defects if they are in circular cylinder shapes. However, the underlying physics for excited modes inside defects is wrongly ignored. In this work, we uncover that additional modes could be excited by theoretical analysis, which is important as it will correct the current common perception. For the case of matched impedance zero-index metamaterials (MIZIM), the additional dipole modes can be excited inside the defects when total transmission occurs. Moreover, we also observe the same results in Dirac-cone-like photonic crystals which have been demonstrated theoretically and experimentally to function as MIZIM. For another case of epsilon-near-zero metamaterials (ENZ), we find that additional higher order modes (e.g., tri-pole) can be excited inside the defects when total transmission happens. Numerical simulations are performed to verify our finding regarding the additional modes.

I n the beginning of the new century, the first metamaterial was achieved to realize the function of negative refractive index 1 . After that, the research on metamaterials 2-7 has made great progress. Recently, attention to zero-index-metamaterials (ZIM) [32][33] has been extensive. For instance, matched impedance zero-index metamaterials (MIZIMs), epsilon-near-zero metamaterials (ENZ), anisotropic ENZ. By utilizing ZIM, some applications and devices with novel functionalities can be realized, such as squeezing wave energy [10][11][12][13][14] , tailoring wave front [15][16][17] , realizing total transmission and reflection in ZIM [18][19][20][21][22] , waveguide bending 23 , enhancing radiation from an embedded source [24][25][26] , controlling energy flux 27 , etc. Several years ago, by putting perfect electric conductor (PEC) or perfect magnetic conductor (PMC) defects in ZIM in a waveguide structure, Hao et al. 18 confirmed that incident electromagnetic wave can undergo total reflection or transmission. Later, Nguyen et al. 19 found similar effects by introducing dielectric defects into MIZIM. However, due to an insufficient expression of the magnetic field, some interesting phenomena and physics are thereby missing. For example, they claimed that for total transmission, only monopole modes exist in the dielectric defects. A lot of publications [20][21][22][28][29][30] (including one from the authors 20 ) later followed this erroneous step, albeit some other intriguing properties found. In this letter, we will give a more comprehensive analysis and show that additional higher modes are excited together with the monopole modes. Our paper corrects some common misunderstanding and shows more colorful physics for ZIM systems.

Results
Now let us start from the schematic plot of a two dimensional (2D) waveguide structure in Fig. 1. Region 0 and region 3 are free space. Region 1 is ZIM with the effective permittivity and permeability e 1 and m 1 . Region 2 consists of N cylindrical defects embedded in region 1. The effective permittivity and permeability of the j-th cylinder are e 2j and m 2j , respectively. Without loss of the generality, we suppose that a transverse magnetic (TM) wave (its magnetic field H is along z direction) is incident from the left port of the waveguide. The outer boundaries of the waveguide are set as PECs. If the incident wave is a transverse electric (TE) wave (with the electric field E polarized along z direction), the outer boundaries of the waveguide should be changed into PMCs for similar results.
For simplicity, we assume the incident magnetic field H I int~ẑ H 0z e i k0x{vt ð Þ , where k 0 is the wave vector in free space with k 0 5 v/c, v is the angular frequency, c is the velocity of light in free space, H 0z is the amplitude of the OPEN SUBJECT AREAS: MICRO-OPTICS SUB-WAVELENGTH OPTICS incident magnetic field. In the following sections, we will omit the time harmonic factor e 2ivt . The electromagnetic (EM) wave in each region satisfies the Maxwell's equations: where the integer m indicates each region and e m is the relative permittivity of each region. The magnetic field in region 0 is a summation of the incident wave and the reflected wave, and is written as, the electric field is, where < is the reflection coefficient. Likewise, we can obtain the magnetic field and the electric field in region 3 as, where = is the transmission coefficient. In region 1, as e 1 almost equals to zero, +|H I 1 must be zero in order to guarantee a finite E 1 . Consequently, the magnetic field in region 1 denoted as H 1 should be a constant. By applying the boundary conditions at the interfaces of x 5 0 and x 5 a, we have, Therefore, <z1~=: In region 2, the magnetic field inside each cylindrical defect follows the Helmholtz equation, The solution can be written as a summation of infinite number of Bessel functions with angular terms. Therefore, the magnetic field in region 2 should be written as, where t jn are the coefficients to be determined for the n-th order Bessel functions.
By applying Dirichlet boundary conditions at the surface of each defect, the magnetic field in region 2 becomes, where J n (x) is the n-th order Bessel function with J n (k 2j R j ) 5 0. a jn and b jn are coefficients of the excited higher order modes, k 2j~k0 ffiffiffiffiffiffiffiffiffiffi e 2j m 2j p is the wave vector in the j-th cylindrical defect, R j is the radius of the j-th cylinder, r j is relative radial coordinate in the j-th cylinder, h j is relative angular coordinate in the j-th cylinder, as mentioned in Ref. 19. We note that the original expression of the magnetic field in each defect is wrongly assumed in Ref. 19 from the two missing terms that typify the additional modes in the defects. However, if J n (k 2j R j ) ? 0, the two additional terms should not be included so that at each circular boundary r j 5 R j , the magnetic field takes the constant value H 1 .
With Eq. (9), the electric field inside each defect could be obtained as follows by recalling Eq. (1), whereĥ j is the azimuthal unit vector for the j-th cylindrical defect. By using the Maxwell-Faraday equation, and after some meticulous calculations, we could obtain the transmission coefficient as 20 , where S 5 a 3 h is the entire area of region 1 and region 2, pR 2 j is the total area of region 2, which consists of N cylinders.
For the system of MIZIM For matched impedance zero-index materials (MIZIM), m 1 %0, Eq. (12) changes into the following formula, which has been derived in Ref. 19 (the Eq. (8) therein). From the Eq. (13), we see that to achieve total transmission (=?1), J 1 (k 2j R j ) must be equal to zero. Therefore, we should select n 5 1 in Eq. (9), the magnetic field in region 2 is then written as, It is not only a function of r j but also a function of h j , a j1 and b j1 are coefficients to be determined. Eq. (14) is the unique solution for total transmission with J 1 (k 2j R j ) 5 0. According to Eq. (14), we note that the magnetic field inside each defect consists of not only monopole modes, but also additional dipole modes.
To solve the coefficients of dipole modes, we suggest an approximate system, i.e., a dielectric cylinder embedded in MIZIM as the background. As the waveguide system supports TM 0 mode, which resembles a plane wave, we suppose the incident wave is a plane wave along x-direction. The magnetic field in MIZIM could be expressed as 31 , where H n (x) is the n-th order Hankel function of the first kind, H 1 takes the same value of the constant magnetic field in the MIZIM area, g n is the scattering coefficients, k 1~k0 ffiffiffiffiffiffiffiffi ffi e 1 m 1 p is the wave vector in MIZIM. When total transmission occurs, the scattering coefficients should be very tinny, that is g n < 0.
The magnetic field in the dielectric cylinder could be written as, At the boundary of the cylinder r 5 R 1 , we should have, We can easily obtain the coefficient of each mode, By combining Eq. (16) and (18) with Eq. (9), we could obtain the relationship between t 1n and a 1n as follows, For MIZIM, as Eq. (14) only has two terms (n 5 0 and n 5 1). We could therefore obtain that, a 0 is the coefficient of the monopole mode, and a 11 is the coefficient of the dipole mode. Due to the symmetry in y-direction, the degenerate state of dipole mode inside the defect (J n (k 2j r j )sin(nh j ) in Eq. (9)) could not be excited out, the coefficient of the degenerate state (b jn ) should be zero. We plot the relationship between the parameter ja 11 j/ja 0 j and different frequencies for three types of MIZIM (with e 1 5 m 1 5 10 23 , 10 24 and 10 25 ) from Eq. (21) and (22), as shown in Fig. 2. We also plot the related numerical result from the waveguide structure (the simulation results in the following sections are all from COMSOL). We find that dipole mode is much more dominative than monopole mode at the resonance frequency that J 1 (k 21 R 1 ) 5 0. While for the frequencies slightly away from the resonance frequency, dipole mode becomes very weak or diminishes. When e 1 5 m 1 tends to zero gradually, the coefficient of dipole mode will decrease accordingly, and the resonance peak becomes narrower. The resonance frequencies from the numerical results have tiny shifts due to the effect of PEC boundaries of the waveguide. When e 1 5 m 1 tends to zero gradually, they will get closer to the analytic resonance frequency.
We will verify the above findings from the numerical simulations. For simplicity, we assume that there is only one cylindrical defect in the region of ZIM. The radius of the cylindrical defect is 0.2 m. Its dielectric constant is e 21 5 4, and its relative permeability is m 21 5 1. In order to achieve total transmission, the term J 1 (k 21 R 1 ) should be equal to zero. As a result, the working frequency is 0.45737 GHz. However, due to the effect of PEC boundaries of the waveguide, this frequency is slightly away from the real resonant frequency, leading to diminishing dipole mode. Therefore, in the simulations, we select the frequency of 0.457 GHz so as to obtain more clear dipole mode. In the region of ZIM, we set a 5 h 5 0.8 m. The magnetic field distributions for the above waveguide system are shown in Fig. 3 for ZIM with different permittivities. In Fig. 3(a), we set e 1 5 m 1 5 0.001 and H 1 5 1, and it seems that only monopole mode is excited from the field pattern. However, the field pattern just takes the real part of the magnetic field in COMSOL. In fact, if we read from the imaginary part, we could find that dipole mode exists. The magnetic field inside the defect is a summation of monopole mode and dipole mode from Eq. (14). The coefficients of monopole mode and dipole mode have a p/2 phase difference from Eq. (21) and Eq. (22). For instance, if we change H 1 5 1 into H 1 5 i, we could find that dipole mode appears from the field pattern in Fig. 3(b). If we set e 1 5 m 1 5 0.0001, the dipole mode will become weaker, as shown in Fig. 3(c). For e 1 5 m 1 5 0.00001, dipole mode almost disappears, see in Fig. 3(d). We can also see this from Fig. 2. Numerically, we find that the coefficients of dipole modes are a 11 5 12.37i, a 11 5 23.42i and a 11 5 20.245i for the cases of Fig. 3 (b), Fig. 3 (c) and Fig. 3 (d), respectively. For all the above cases, the coefficients of the monopole modes are 1/J 0 (k 21 R 1 ) 5 22.483 and the coefficients of the degenerate state of dipole modes are zero (b 11 5 0). For e 1 5 m 1 5 0.001, dipole mode is more dominative than monopole mode. While for e 1 5 m 1 5 0.0001, they are comparable to each other. For e 1 5 m 1 5 0.00001, the monopole mode is more dominative than dipole mode. It seems that the dipole mode is diminishing when e 1 5 m 1 tends to zero gradually. It is not easy for us to choose the required resonant frequency as the resonance goes extremely narrow. For a particular value of near zero permittivity, we can in principle find a frequency  near the resonant one where dipole mode is much more dominative than monopole mode.
To demonstrate the hybridization of monopole mode and dipole mode more clearly, we plot the real part of magnetic field distribution inside the defect from x 5 20.2 m to x 5 0.2 m (at y 5 0) with different values of H 1 for e 1 5 m 1 5 0.001, as shown in Fig. 4. For H 1 5 1, only the information of monopole mode is observed from the black curve (it is an even function of x). While for H 1 5 i, the information of dipole mode could be observed from the red curve (it is an odd function of x). For H 1 5 0.707 1 0.707i, both information of monopole mode and dipole mode could be observed from the blue curve. To make it more straightforward, we also plot the amplitude of the magnetic field (see the green curve), which is independent of H 1 and has two symmetric peaks at the positions of x 5 0.1 m and x 5 20.1 m because of the existence of dipole mode. If there is only dipole mode inside the defect, the amplitude should be zero at the position of x 5 0. However, the amplitude there is a value of about 2.5, which is equal to the amplitude of the monopole mode at the position of x 5 0 (see also the black curve). Therefore, both monopole mode and dipole mode exist inside the defect.
As it is known, Dirac-cone-like photonic crystals 16 can be regarded as MIZIM near Dirac point frequency. It should be possible to pro-duce the above similar effect if we replace MIZIM with such photonic crystals. It is noticed that the incident wave is now a transverse electric (TE) wave, and the outer boundaries of the waveguide structure are PMCs. We will show that if a cylindrical defect is introduced in MIZIM or Dirac-cone-like photonic crystals, at the condition of J 1 (k 21 R 1 ) 5 0 when total transmission occurs, both dipole mode and monopole mode exist inside the defect in the waveguide system (k 21~k0 ffiffiffiffiffiffiffiffiffiffiffi ffi e 21 m 21 p is the wave vector of light in the defect, e 21 and m 21 are the permittivity and permeability of the defect respectively, R 1 is the radius of the defect). Following Ref. 16, the Dirac-cone-like photonic crystals consist of cylindrical alumina rods arranged in a square lattice. The radii of the rods are 3.75 mm with a dielectric constant 8.8. The lattice constant is 17 mm. The Dirac point frequency f is about 10.3 GHz. We set a 5 h 5 0.187 m for the region of photonic crystals in the waveguide and insert a cylindrical defect with a radius of R 1 5 0.0284 m in the center of the region (the radius should be large enough to visualize the above effect). The dielectric constant of the defect is e 21 5 0.3905 and its relative permeability is m 21 5 1, to satisfy J 1 (k 21 R 1 ) 5 0 at the Dirac point frequency. In Fig. 5 (a), we plot the electric field for the system of Dirac-cone-like photonic crystals and choose a suitable phase of the incident plane wave (E 1 5 i) so that only dipole mode is demonstrated in the cylindrical defect. Likewise, we choose another phase of the incident plane wave (E 1 5 1), and plot the electric field in Fig. 5 (b), where only monopole mode is shown in the defect. For comparison, we replace the Diraccone-like photonic crystals with MIZIM, and plot corresponding electric field in Fig. 5 (c) and (d). Fig. 5 (c) shows a consistent dipole mode with Fig. 5(a), while Fig. 5(d) gives out monopole mode like that in Fig. 5(b). Therefore we cannot neglect the existence of dipole mode like Ref. 19. Sometime, it is more dominative than the monopole term near the resonance frequency, even for such a realistic photonic crystal system.

For the system of ENZ
After discussing the MIZIM case, we come to the ENZ case. Let us return to Eq. (12). In order to obtain total transmission, the following term must be zero, For simplicity, we suppose that there is only one cylindrical defect in the ENZ area. After some calculations, we shall have, If the n-th order Bessel function J n (k 21 R 1 ) is zero, the magnetic field in the defect could be written as, Therefore, there is not only monopole mode excited inside the defect, but also some other additional higher order mode emerging as well, if both the conditions of Eq. (24) and J n (k 21 R 1 ) 5 0 are satisfied. For example, we choose J 3 (k 21 R 1 ) 5 0 and tune the configuration and material parameters to satisfy Eq. (24). Following similar calculations to Eq. (19) and (20), we could get the coefficient of tripole mode a 13 5 i 3 2J 3 (k 1 R 1 )/J 3 (k 21 R 1 ) and the coefficient of monopole mode a 0 5 1/J 0 (k 21 R 1 ). Likewise, we find the relationship between the parameter ja 13 j/ja 0 j and different frequencies for three types of ENZ (with e 1 5 10 23 , 10 24 and 10 25 ) both theoretically and numerically, as shown in Fig. 6. From the analytical results (solid curves), we find that tri-pole mode is much more dominative than monopole mode at the resonant frequency where J 3 (k 21 R 1 ) 5 0. For other frequencies slightly deviating from the resonant one, tri-pole mode is disappearing. When e 1 tends to zero gradually, the coefficient of tri-pole mode will decrease accordingly and the resonance peak will become narrower. However, the resonance frequencies from numerical results have a tiny shift because of the effect of outer PEC boundaries of the waveguide. When e 1 tends to zero gradually, the resonance frequency will approach the analytic resonance frequency.
The finding will be confirmed again from numerical simulations. Suppose that there is one cylindrical defect inside the ENZ area. The radius of the defect is 0.01 m, its dielectric constant is e 21 5 16 and its relative permeability is m 21 5 1. In order to make J 3 (k 21 R 1 ) 5 0, where k 21 R 1 is the first root of the third order of Bessel function, the working frequency should be about 7.615 GHz. In addition, to satisfy Eq. (24), we set the effective permittivity and permeability of the ENZ as e 1 5 0.001 and m 1 5 0.6 respectively and set a 5 h 5 0.021 m for the ENZ area. The magnetic field distribution for the above system is shown in Fig. 7. The tri-pole mode is demonstrated inside the defect when H 1 5 i, as shown in Fig. 7(a). Likewise, we can also observe the monopole mode by changing H 1 into 1, as shown in Fig. 7(b). We numerically find that the coefficient of tri-pole mode is a 13 5 21.67i, b 13 5 0 and the coefficient of the monopole is 1/J 0 (k 21 R 1 ) 5 4.17 (the tripole mode here is obvious but not dominative, to get a dominative tri-pole mode, the working frequency should be shift to about   7.61 GHz, as already shown in Fig. 6). Hence, the simulation results prove our finding, and the magnetic field is a summation of monopole mode and tri-pole mode, which should be written as,