Dynamic study and applications of metamaterial systems

Abstract

We investigate the dynamic characteristics of metamaterial systems, such as the temporal coherence gain of the superlens, the causality limitation on the ideal cloaking systems, the relaxation process and essential elements in the dispersive cloaking systems, and the extending of the working frequency range of cloaking systems. The key point of our study is the physical dispersive properties of metamaterials, which are well-known to be intrinsically strongly dispersive. With physical dispersion, new physical pictures can be obtained for the waves propagating inside metamaterial, such as the “group retarded time” for waves inside the superlens and cloak, the causality limitation on real metamaterial systems, and the essential elements for design optimization. Therefore, we believe the dynamic study of metamaterials will be an important direction for further research. All theoretical derivations and conclusions are demonstrated by powerful finite-difference time-domain simulations.

Graphical abstract

Keywords

metamaterial / dynamic / dispersion

1 Introduction

Materials whose permittivity \epsilon and permeability μ are simultaneously negative are found to possess a negative refractive index n with many unusual properties [1]. Negative-n metallic resonating composites and two-dimensional (2-D) isotropic negative-n material have been constructed [2, 3], and negative light refraction was observed [4]. The unconventional properties of such materials, such as the evanescent wave being able to be amplified by negative-n so that the sub-wavelength resolution can be achieved [11, 12], have drawn an increasing amount of attention in both science and engineering [7]. After negative-n material, more such unconventional materials are found, so that a new concept “metamaterial” is generated, which is an effective medium both for the permittivity {\epsilon }_{\mathrm{per}} and for the effective permeability {\epsilon }_{\mathrm{eff}} over a certain finite frequency band. Such physical media are composed of distinct elements (photonic atoms) which are generally constructed by metallic material, and their size and spacing are on a scale much smaller than the wavelengths in the frequency range of interest. Thus, the effective composite media can be considered homogeneous at the wavelengths under consideration. Because of their abnormal properties which can even go beyond the traditional physical limit and totally new optical phenomena, metamaterial has become one of the hottest topics in modern photonics. However, many of the current studies on metamaterials are mainly concentrated on the stationary state, in which only the single-frequency properties are considered in the study. It is understandable since the stationary state is most important because it represents the stable working state of metamaterial, and the studies of stationary state are much easier than others, i.e., the dynamic studies. However, without dynamic studies, many basic questions of stationary state cannot be answered, such as, “How can the field get to its stable state?”, “Is there any strong scattering or oscillation in the process?”, “How long is the relaxation process?” and “What is the system response if the incidence is a pulse?”. Furthermore, some physical topics that are intrinsically dynamic, for example, the temporal coherence change in propagation, the tuning of cavity frequency, and the nonlinear response, etc.. The essential point of the difference between the dynamic study and the stationary study is whether the dispersion of metamaterial is neglected. For common dielectric material, if the dispersion is very weak in the frequency range of interest, the stationary study can give us a pretty clear physical picture to understand the optical properties. However, we know that metamaterials are usually highly dispersive in the window of the work frequency, i.e., for the thin wire medium the permittivity can be described by the effective dielectric function {\epsilon }_{\mathrm{eff}}=1-{\omega }_p^2/{\omega }^2 with the plasma frequency {\omega }_p related to the geometry of the wire array with typical value {\omega }_p\sim 8.2\mathrm{GHz} and for the SRR medium the effective frequency-dependent permeability, usually having the form {\mu }_{\mathrm{eff}}=1-F{\omega }_0^2/({\omega }^2-{\omega }_0^2-\mathrm{i}\omega \Gamma ) [8]. The deeper physical reason is that the strong interaction between electromagnetic field and “photonic atoms” is the origin of almost all abnormal optical properties of metamaterial, and it is the origin of the strong dispersion too. In other words, strong dispersion and the abnormal optical property are twins with the same origin. Therefore, it is not a surprise that the dispersion is generally needed to explain those abnormal optical properties, such as the negative refraction of negative-n material. More seriously, if the specific strong dispersion of a certain metamaterial is neglected, the violation of the basic physical laws is pointed out, i.e., it can be deduced that there exists causality violation, such as the superluminal group velocity in cloaking systems [9] and hyper-lens systems [10]. Hence, the study with dispersion of metamaterial is not only needed, but also required. Although for a linear system we can obtain the dynamic picture theoretically if we know the response of the system at all frequency range for the strongly dispersive metamaterial, is so hard that the dynamic response needs to be studied specifically. In these dynamic studies, the numerical simulation is an important tool to obtain direct observation of the details. The well-developed finite-difference time-domain (FDTD) method is a good candidate for the dynamic simulation of metamaterials.

In this paper, we will emphasize the dynamic picture of metamaterial systems, such as the superlens and the electromagnetic cloak. As we have discussed, the dispersion of metamaterial is the key to the whole study. From dynamic studies, we will not only reveal the deeper physical pictures of metamaterial systems, but also discuss the physical limitation in these systems. The paper is organized as follows.

In Section 2, the image field of the negative-index superlens with a quasi-monochromatic random source is discussed, and dramatic temporal-coherence gain of the image in the numerical simulation is observed, even if there is almost no reflection and no frequency filtering effects. From the new physical picture, a theory is constructed to obtain the image field and demonstrate that the temporal coherence gain is from different “group” retarded time of different optical paths. Our theory agrees well with the numerical simulation and strict Green’s function method. This study should have important consequences in the coherence studies in the related systems and the design of novel devices.

In Section 3, the dynamic processes of the dispersive cloak by finite-difference time-domain numerical simulation are carried out. It is found that there is a strong scattering process before the stable state is achieved and its time length can be tuned by the dispersive strength. Poynting-vector directions show that the stable cloaking state is constructed locally while an intensity front sweeps through the cloak. Deeper studies demonstrate that the group velocity tangent component V_{g\theta } is the dominant factor in the process. This study is helpful not only for clear physical pictures but also for designing better cloaks to defend passive radars.

In Section 4, the limitation of the electromagnetic cloak with dispersive material is investigated based on causality. The results show that perfect invisibility cannot be achieved because of the dilemma that either the group velocity V_g, diverges or a strong absorption is imposed on the cloaking material. It is an intrinsic conflict which originates from the demand of causality. However, the total cross section can really be reduced through the approach of coordinate transformation. A simulation of FDTD method is performed to validate the analysis.

In Section 5, we investigate the optical properties of the interface between the anisotropic metamaterial with hyperbolic dispersion and the isotropic dielectric. With material dispersion, a comprehensive theory is constructed, and the hyperlens effect that the evanescent wave can be converted into the radiative wave is confirmed. At the inverse process of the hyperlens, we find a novel mechanism to compress and stop (slow) light at wide frequency range, which can be used as a removable memory or a light trap.

In Section 6, based on the transformation optics, a general method of light-behavior remote control is proposed. From this method, the important coefficients of a cavity, i.e., the quality factor Q and the resonant frequency {\omega }_0 can be tuned in a wide range by a transformation optical device in distance, so that the light behavior can be remotely controlled. To confirm this original idea, three schemes, such as the remote modification of output energy current from an absorptive cavity, the remote control of lasing behaviors, and the remote tuning of the resonant frequency or photonic band-gap, are presented and confirmed by our numerical simulations based on finite-difference time-domain and finite-element methods. With some special advantages, e.g., without physical change or damage of original devices, large tuning range, and easily hidden controller, this method can be widely used in optical/photonic or electromagnetic designs in the future.

2 The temporal coherence gain of the negative-index superlens image

Veselago predicted that the negative-index material (NIM) has some unusual properties, for instance, the flat slab of the NIM could function as a lens for electromagnetic (EM) waves [1]. This research direction was further pushed by works of Pendry and others [4-6, 11-28] who showed that the lens with such NIM (i.e., \epsilon =\mu =-1+\delta) could be a superlens whose image resolution can go beyond the usual diffraction limit. After that, several beyond-limit properties of NIM systems are found, such as, the sub-wavelength cavity [29] and the waveguide [30, 31]. Some of the theoretical results are confirmed by experiments [4-6, 16, 17]. Also, these beyond-limit properties give us new physical pictures and opportunities to design devices. Recently, new numerical [18-20] and theoretical Green’s function [21, 22] methods were used to understand the phenomena in such systems. But so far almost all studies are done with the strictly single-frequency sources, so that the coherent properties of EM waves (or photons) in the NIM systems have not been studied to the best of our knowledge. Even more seriously, there is no theory for the propagation of coherent functions in NIM systems. The importance of coherence research cannot be over-estimated since coherence is essential in wave interference, imaging, signal processing and the telecommunication [32-34]. Can we find a new frontier to go beyond at the coherent properties in NIM systems? If so, can we develop a simple theoretical method to deal with the image coherence of the superlens?

In this section, the FDTD method is used in the 2-D numerical experiments to study the temporal coherence of the superlens image with random quasi-monochromatic sources. We observe the dramatic temporal-coherence gain of the superlens image even if the reflection and frequency-filtering effects are very weak. Basing on the new physical picture of the signal (the fluctuation of random source) propagation in NIM, we construct a theory to obtain the image field and derive the equation of the temporal-coherence relation between the source and its image. The new mechanism of the temporal-coherent gain can be explained by the key idea that the signals on different paths have different “group” retarded time. Our theory agrees completely with numerical results and the strict Green’s function results.

The setup of the 2-D system is shown in Fig. 1. The thickness of the infinitely-long NIM slab is d. To realize the negative \epsilon and negative μ, the electric polarization density P and the magnetic moment density M are phenomenologically introduced in FDTD simulation [35]. The effective permittivity and permeability of the NIM are {\epsilon }_r(\omega )={\mu }_r(\omega )=1+{\omega }_P^2/({\omega }_a^2-{\omega }^2-\mathrm{i}\gamma ). In our model, {\omega }_a=1.884\times {10}^{13}/\mathrm{s}, \gamma ={\omega }_a/100, {\omega }_P=10\times {\omega }_a. The quasi-monochromatic field is expressed as E(x,t)=U(x,t)\exp (-\mathrm{i}{\omega }_{0}t), where U(x,t) is a slowly-varying random function, {\omega }_0\text{π}=/20{\delta }_t is the central frequency of our random sources and {\delta }_t=1.18\times {10}^{-15}s is the smallest time-step in FDTD simulation. At {\omega }_0, we have {\epsilon }_r={\mu }_r=-1.00-0.0029\mathrm{i}. Here, we emphasize that in our FDTD simulation the smallest space-step {\delta }_x={\lambda }_0/20({\lambda }_0=2\text{π}c/{\omega }_0) and the distance (d/2={\lambda }_0) of the source from the lens are too large to excite strong evanescent modes of NIM [18-20, 23]. Actually the evanescent field in our simulation can be neglected compared with radiating field, and what we are studying is the property dominated by the radiating field.

The random source is composed of the randomly generated plane-wave pulses, with the average pulse length t_p and the random starting phase and starting time. In the simulation, we record the field of the source and the image for a duration of 4\times {10}^5{\delta }_t to obtain the data for analysis. For convenience, we define E(\omega )={lim}_{T\to \mathrm{\infty }}{\int }_{-T}^{T}E(t)\exp (-\mathrm{i}\omega t) as the field spectrum (FS).

Unusual phenomena.-At first, the FS width of the random source is a bit too large (\mathrm{\Delta }{\omega }_s\simeq {\omega }_0/20). When we observe the image temporal-coherence gain, we also find that the FS width of the image is sharper than the source (). It is obvious that there are frequency-filtering effects because of the NIM dispersion, such as frequency-dependent interface reflection and focal length. After increasing the pulse-length t_p of the source, we reduce the source FS width to \mathrm{\Delta }{\omega }_s\simeq {\omega }_0/100, then the reflection and focal-length difference are very small [36]. With such source, the FS widths of source and image are almost the same \mathrm{\Delta }{\omega }_i\simeq \mathrm{\Delta }{\omega }_s, as shown in Fig. 2(a). The difference between the two widths is , which is our criterion of the quasi-monochromatic source. Even so, the dramatic gain of temporal coherence is still observed. In Fig. 2(b), the source field (up) and the image field (down) vs. time of FDTD simulation are compared. Their profiles are generally similar, but the image profile is much smoother. The normalized temporal-coherence function g^{(1)}(\tau )=⟨E^{\ast }(t)E(t+\tau )⟩/⟨E^{\ast }(t)E(t)⟩(⟨⟩ means the ensemble average) of the source (black) and the image (red) from FDTD simulation are shown in Fig. 4. The temporal coherence of the image field is obviously better than the source. From g^{(1)}, the image coherent time is obtained T_i^{\mathrm{co}}=\int g_i^{(1)}(\tau )d\tau =1268{\delta }_t, which is about 50\% longer than the source coherent time T_s^{\mathrm{co}}=860{\delta }_t.

Although the gain of the spatial coherence only by propagation is well-known [33, 34], the dramatic gain of temporal coherence is generally from the high-Q cavities, contrary to our case, which have strong filtering effects. To reveal the new mechanism of the temporal coherence gain in NIM systems, we also have done more numerical experiments in which only the ray near a certain incident angle (shown in Fig. 1), such as only paraxial rays (\theta \simeq 0), can pass through the superlens. Then the image field profile vs. time looks much like the source field and has no gain of coherence anymore. Therefore, the gain of temporal coherence of the superlens image is not from one ray with certain incident angle, but probably from the interference between the rays with different incident angles. Then, what is different between the rays with different incident angles? After carefully checking the field profiles of different-incident-angle cases, we find that the profiles have different retarded time. The larger the incident angle the longer the retarded time.

Physical pictures. -To deeply understand the new mechanism of coherence gain and construct our theory, we  need  to  make  two physical pictures clear. The first one is about the optical path length (OPL) \int nds which determines the wave phase and the refracted “paths” of rays in Fig. 1 according to Fermat’s principle (or Snell’s law). Based on ray optics, the superlens and traditional lenses have the same focusing mechanism, in which all focusing rays have the same OPL {\int }_{\mathrm{paths}}n\mathrm{d}s=const({\int }_{\mathrm{paths}}n\mathrm{d}s=0 for superlens) from source to image [1]. However, this picture is so well-known that it suppresses the other important picture. Because the temporal-coherence information is in the fluctuation signals of the random field, the signal propagating picture should be essential to our study. The optical signals propagate in the group velocity v_g which is always positive. Obviously, if the path (in Fig. 1) is longer (larger incident angle), the signal needs a longer propagating time, which is called group retarded time (GRT) in this section. Inside the NIM, the GRT of a path should be d/[\cos (\theta )v_g], (which is confirmed by our numerical experiments,) where θ is the incident angle and v_g=c/3.04 is the group velocity of NIM around {\omega }_0 [37]. The total GRT from source to image is {\tau }_r={\tau }_0/\cos \theta where the {\tau }_0=d/c+d/v_g is the GRT of the paraxial ray. Now, the new propagating picture for a signal through superlens is that a signal, generated at t_s from the source, will propagate on all focusing paths and arrive at the image position at very different time t_s+{\tau }_0/\cos \theta from different paths (this is schematically shown in Fig. 1). This picture is totally different from traditional lenses, whose images do not have obvious temporal-coherence gain because their focusing rays have the same OPL and similar GRT.

Our theory. -Based on these analysis, we suppose that the superlens image field of the random quasi-monochromatic source is the sum of all signals from different paths with different GRT. This is the key point of our theory, and then the image field can be obtained:

\begin{array}{}\text{(1)}& \begin{array}{cc}{E}_i(t)\hfill & =\frac{1}{U_0}{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}{\displaystyle \sum _{\mathrm{paths}}}{U}_s(t-{\tau }_r)\hfill \\ \hfill & =\frac{1}{U_0}{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}{\displaystyle {\int }_{-\frac{\text{π}}{2}}^{\frac{\text{π}}{2}}}{U}_s(t-\frac{{\tau }_0}{\cos \theta })d\theta \hfill \end{array}\end{array}

where U_s(t) is the slowly-varying profile function of the source and U_0 is the normalization factor. In Fig. 2(c) (up), we show the result of the image field based on Eq. (1), and we can see it is in good agreement with the FDTD result in Fig. 2(b) (down). To show the interference effect of different paths, we assume there are only two paths (such as A and B in Fig. 1). Based on Eq. (1) the image field is E_i={\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}[U_s(t-{\tau }_r^A)+U_s(t-{\tau }_r^B)], then the temporal coherence of image is G(\tau )=⟨E_i^{\ast }(t)E_i(t+\tau )⟩=⟨U_s^{\ast }(t-{\tau }_r^A)U_s(t-{\tau }_r^A+\tau )+U_s^{\ast }(t-{\tau }_r^B)U_s(t-{\tau }_r^B+\tau )+U_s^{\ast }(t-{\tau }_r^A)U_s(t-{\tau }_r^B+\tau )+U_s^{\ast }(t-{\tau }_r^B)U_s(t-{\tau }_r^A+\tau )⟩. The first two terms are the same as the source field (just a time-shift) so they do not contribute to the coherence gain. The last two terms are from interference between two paths. The third (or the fourth) term could be very large at the condition \tau \simeq {\tau }_r^B-{\tau }_r^A (or {\tau }_r^A-{\tau }_r^B). This condition can always be satisfied between any two paths since τ is a continuous variable. So the interfering terms between the paths are responsible for the image temporal-coherence gain.

From Eq. (1), after the variable transformation t_s=t-{\tau }_0/\cos \theta and some algebra, the relation of the temporal coherence between the image and the source can be obtained:

\begin{array}{}\text{(2)}& \begin{array}{cc}{G}_i(\tau )\hfill & =⟨E_i^{\ast }(t)E_i(t+\tau )⟩\hfill \\ \hfill & =\frac{1}{U_0^2}{\displaystyle {\int }_{-\mathrm{\infty }}^{-{\tau }_0}}d{t}_1{\displaystyle {\int }_{-\mathrm{\infty }}^{-{\tau }_0+\tau }}\text{d}{t}_2{h}_i^{\ast }(t_1)\cdot h_i(t_2+\tau )G_s(t_2-t_1)\hfill \end{array}\end{array}

where h_i(t)={({\tau }_0/t)}^2/\sqrt{1-{({\tau }_0/t)}^2} is the response function of different incident angles, and G_s(t_2-t_1)=⟨E_s^{\ast }(t_1)E_s(t_2)⟩ is the temporal-coherence function of the source. Equation (2) can explain the temporal-coherence gain of the image too. Even if the source field is totally temporal incoherent G_s(t_2-t_1)\propto \delta (t_2-t_1), based on Eq. (2), we can find that G_i(\tau ) is not a δ-function anymore, so the image is partially temporal coherent. The product of h_i^{\ast }(t_1)h_i(t_2+\tau ) includes the interference between paths. According to our theory, we calculate the coherent T vs. the superlens length L (Fig. 3 blue) which agrees with our FDTD result (Fig. 3 red) pretty well. (We will discuss the deviation later.)

To further confirm our theory and FDTD results, the strict Green’s function method [21, 22] is engaged to check our results. We only include the radiating field (no evanescent wave) in Green’s function. The strict image field vs. time is shown in Fig. 2(c) (down), and the image temporal-coherence function g^{(1)} vs. time is shown in Fig. 4 (blue). In Fig. 4, we can see that the FDTD result (red) is almost exactly the same as the strict Green’s function method (green). However, our theory (blue) deviates from the strict result at very large \tau >3000{\delta }_t which corresponds to very long path (or very large incident angle). This is understandable since in our theory we neglect the dispersion of NIM totally and use only v_g({\omega }_0). For the very-large-angle rays a small index difference (from the dispersion of NIM) can cause large focal-length difference. Hence, the deviation is from the focus-filtering effect. When we reduce the FS width of source to an even smaller value (i.e., \mathrm{\Delta }{\omega }_s={\omega }_0/500), the deviation of our theory is smaller.

Although our theory is merely a good approximation generally, owing to the picture simplicity and clarity, the theory can help us to study more complex systems qualitatively and quantitatively. The finitely-long 2-D superlens is a good example which is hard to deal with by Green’s function method. In Fig. 4, we plot the coherent time T_i^{\mathrm{co}} vs. superlens length L of the FDTD simulation (black) and of our theory (red), respectively. They coincide with each other pretty well. (The deviation reason has been discussed.) The increase of the T_i^{\mathrm{co}} with the increase of L can be explained simply according the our theory. Since the image field is E_i(t)=1/U_0{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}{\int }_{{\theta }_{\min }}^{{\theta }_{\max }}{U}_s\left(t-{\tau }_0/\cos \theta \right)d\theta, the large-angle paths (\theta >{\theta }_{\max } and ) and their contribution to the temporal-coherence gain are missed in the short superlens.

Obviously, Eq. (1) is suitable not only for random quasi-monochromatic source, but also for all quasi-monochromatic fields, such as the slowly-varying Gaussian pulses and slowly switching-on process mentioned in Refs. [21, 22]. Our theory can be easily extended to 3-D systems too. Also, owing to the fact that what we have found is from the radiating field, the temporal-coherence gain is not the near-field property. Actually, the new mechanism of the temporal-coherence gain is not limited to the n\simeq -1 superlens, but also applicable to other superlenses, like the photonic crystal superlens in Refs. [5, 6, 16, 28]. However, the specialities of n\simeq -1 superlens, such as almost no frequency-filtering (no frequency loss) and no reflection (no energy loss), can be used to design novel optical/photonic coherence-gain devices.

In summary, for the first time we have numerically and theoretically studied the temporal coherence of the superlens image with the quasi-monochromatic source. Numerically, we have observed that the temporal coherence of the image can be improved considerably even almost without reflection and filtering effects. Basing on the new physical picture, we construct a theory to calculate the image field and temporal-coherence function, which excellently agree with the FDTD results and strict Green’s function results. The mechanism of the temporal coherence gain is theoretically explained by the different GRT of different paths. Although the evanescent wave is very weak in this study, the coherence of evanescent wave in NIM systems is an interesting topic which will be discussed elsewhere [38]. Other related topics, such as the spatial coherence which is essential to the image quality of the superlens, can also be studied through the similar methods. Although our study is within the confinement of classic optics, similar investigation can be extended to the quantum optics [33, 34], and interesting results can be expected. Obviously, the temporal-coherence gain of superlens is another evidence that the NIM phenomena are consistent with the causality [16]. We suppose that the temporal-coherence gain phenomena could be observed in micro-wave experiments [4-6]. Therefore, this study should have important consequences in the future studies of coherence in NIM systems. The no-reflection and no-frequency-filtering coherence gain of the superlens has some potential applications in imaging, coherent optical communication, and signal processing.

3 The physical picture and essential elements of the dynamic process for dispersive cloaking structures

Recently, a theory [39, 40] has been developed based on the geometry transformation to realize a cloaking structure (CS), in which objects become invisible from outside. Then a 2-D cylindrical CS [88] and a nonmagnetic optical CS [42, 43] are designed. More surprisingly, the experiment [44] demonstrates that such a 2-D CS really works with a “reduced” design made of split-ring resonators. These pioneers’ works are really attractive and open a new window to realize the dream of human invisibility. However, so far almost all theoretical [45-47, 39, 40, 42, 43, 88] studies of the CS are done in the frequency domain and the geometry transformation idea is supposed to work only for a single frequency, so that the effects of the dispersion have not been intensively studied. As pointed out in Ref. [40] and the quantitative study in our recent work [9], the dispersion is required for the cloaking material to avoid the divergent group velocity. For the dispersive CS, new topics, including the dynamic process, can be introduced. Dynamic study is essential to the cloaking study since without it we cannot answer questions, such as “How can the field get to its stable state?”, “Is there any strong scattering or oscillation in the process?”, and “How long is the process?”, etc.. More importantly, because the real radars are generally pulsive ones, the dynamic process is critical for the cloaking effect around the goal frequency. Therefore, the dynamic study not only gives us a whole physical picture of the cloaking, but also helps us to design more effective cloaks.

In this section, the dynamic process of the electromagnetic (EM) CS is investigated by FDTD numerical experiments. In our simulation, the Lorentzian dispersion relations are introduced into the permittivity and the permeability models, then the real dynamic process can be simulated [48-50]. Based on numerical simulation, we can follow the details of the dynamic process, such as the time-dependent scattered field, the building-up process of the cloaking effect, and the final stable cloaking state. By tuning the dispersion parameters and observing their effects on the dynamic process and the scattered field, we can find the essential elements that dominate the process. Theoretical analysis of these essential elements can help us to have a deeper physical picture beyond the phenomena and to design more effective cloaks.

The setup of the system is shown in Fig. 5(a), similar to the one in Ref. [88]. R₁ and R_2=2R_1 are the inner and the outer cylindrical radii of the CS, respectively. A perfect electric conductor (PEC) shell is pressed against the inner surface of the CS. The CS is surrounded by the free space with {\epsilon }_0={\mu }_0=1. From the left side, an incident plane wave with working frequency {\omega }_0 is scattered by the CS, and the total field and the scattered field can be recorded inside and outside B₁ respectively by the numerical technique [51]. Therefore, the scattering cross-section σ can be calculated easily. Our study is focused on the E-polarized modes, for which only the permittivity and the permeability components {\epsilon }_z, {\mu }_r, and {\mu }_{\theta } are needed to be considered. (For H-polarized modes, considering the corresponding components {\mu }_z, {\epsilon }_r, and {\epsilon }_{\theta }, we can obtain the same numerical results in the dynamic process.) All of them are supposed to have the form 1+F_j(r)\times f_j(\omega ), where subscript j could be z, r, and θ for {\epsilon }_z, {\mu }_r, and {\mu }_{\theta }, respectively. The filling factors F_j(r) are only r-dependent, {\omega }_p is the plasma frequency which is set to be a constant {\omega }_p=10{\omega }_0, and f_j(\omega )={\omega }_p^2/({\omega }_{aj}^2-{\omega }^2-\mathrm{i}\omega \gamma ) are the Lorentzian dispersive functions, where γ is the “resonance width” or called “dissipation factor”, and {\omega }_{aj} are the resonant frequency of “atoms” (resonant units) in metamaterials.

For the study of the dispersive CS, we suppose that the real parts of the {\epsilon }_z, {\mu }_r, and {\mu }_{\theta } always satisfy the geometry transformation of Ref. [88] at {\omega }_0: \mathrm{Re}[{\mu }_r(r,{\omega }_0)]=(r-R_1)/r,\mathrm{Re}[{\mu }_{\theta }(r,{\omega }_0)]=r/(r-R_1), and \mathrm{Re}[{\epsilon }_z(r,{\omega }_0)]=R_2^2(r-R_1)/[{(R_2-R_1)}^{2}r]. Then the  filling  factors   F_j(r)  at  different r can be obtained: F_r(r)=\{\mathrm{Re}[{\mu }_r(r,{\omega }_0)]-1\}/\mathrm{Re}[f_r({\omega }_0)],F_{\theta }(r)=\{\mathrm{Re}[{\mu }_{\theta }(r,{\omega }_0)]-1\}/\mathrm{Re}[f_{\theta }({\omega }_0)], and F_z(r)=\{\mathrm{Re}[{\epsilon }_z(r,{\omega }_0)]-1\}/\mathrm{Re}[f_z({\omega }_0)].

To investigate the dispersive effect on the dynamic process, we tune the dispersion parameters {\omega }_{aj} in our numerical experiments. We use the working frequency {\omega }_0 as the frequency unit since it is the same in all cases in this section so the ratio A_j={\omega }_{aj}/{\omega }_0represents {\omega }_{aj}. Obviously, for the Lorentzian dispersive relation, the dispersion is stronger when {\omega }_0 and {\omega }_{aj} are closer to each other, (The working frequency is near the resonant frequency) or in other words, when A_j approaches one. Since there are singular values of the real parts of \epsilon and μ, we have done some approximations in our numerical simulation [52], for instance, we set the maximum and the minimum for \epsilon and μ. Although such approximations will affect the cloaking effect of the stable state [46], we find that the influence of these approximations on the Dynamic process is very small and can be neglected.

First, we show an example of evolving electronic field during the Dynamic process in Fig. 5 with concrete parameters of A_r, A_{\theta }, A_z, and γ. In Fig. 5(a), the plane wave arrives at the left side of the CS, and is ready to enter the CS. From Fig. 5(b)–(e), the cloaking effect is built up step by step. Finally, the field gets to the stable state shown in Fig. 5(f). Because of the dispersion, there is an obvious time delay in the cloaking effect and a strong scattered field is observed.

We introduce a time-dependent scattering cross-section \sigma (t) to quantitatively study the Dynamic process, which is defined as:

\begin{array}{}\text{(3)}& \sigma (t)={\overline{J}}_{\mathrm{scat}}(t)/{\overline{S}}_{\mathrm{inct}}\end{array}

where t=n\times T,n=0,1,2,\mathrm{\dots }, T is the period of the incident wave, {\overline{J}}_{\mathrm{scat}}(t) is the one-period-average energy flow of scattered field, and {\overline{S}}_{\mathrm{inct}} is the average energy flow density of the incident field. To observe the dispersive effect on \sigma (t) during the Dynamic process, at the first step, we keep A_z and γ constant and change A_r and A_{\theta }, the results are shown in Fig. 6(a). From the σ versus t curves, we can find the general properties of the dynamic process. First, there is strong scattering in the dynamic process. At the beginning, σ increases rapidly when the wave is getting to the CS, then reaches its maximum (at about ninth period). After that, σ starts to decay until it gets to the stable value (of the stable cloaking state). Second, unlike other systems, there is no oscillation in the process. This property will be discussed later. Third, the time length of dynamic process, called “relaxation time” generally, can be tuned by the dispersion. From Fig. 6(a), we can see that the main dispersive effect is on the decaying process. From case 1 to case 5, A_r and A_{\theta } are getting closer to one, so that the dispersion is stronger. We find that the stronger the dispersion, the longer the relaxation time is. To compare with the cloaking cases, we also show the \sigma (t) of the naked PEC shell in the case 6. From the definition of σ, we know that the area covered by these curves in Fig. 6(a) is proportional to the total scattered energy in the dynamic process. Therefore, the CS with the weaker dispersion will scatter less field (better cloaking effect) in the dynamic process. However, such a generalization is still not enough for us to get a clear physical picture to understand the cloaking dynamic process.

Next, we check whether the absorption of the CS is important in the process. The absorption is determined by the imaginary parts of \epsilon and μ. To study this effect, we hold A_r, A_{\theta }, and A_z constant but modify the dissipation factor γ. We modify the filling factors F_j simultaneously, so that the real parts of \epsilon and μ are kept unchanged at {\omega }_0. In such a way, we can keep the dispersion strength almost unchanged, but with the imaginary parts of \epsilon and μ changed. The results in Fig. 6(b) show that the stronger absorption leads only to larger stable value of σ, leaving the relaxation time nearly unchanged. Thus, we can exclude the absorption from the relevant parameter list, since it only influences the \sigma (t) of the stable state considerably.

To obtain a deeper insight of the dynamic process, we need to study it more carefully. From Fig. 5(b)–(e), we can see that the “field intensity” (shown by different colors in the figures) propagates slower inside the CS than that in the outside vacuum. When the inside field intensity “catches up” with the outside one [in Fig. 5(f)], the field in the CS gets to the stable state and the cloaking effect is built up. In fact, this catching-up process of the field intensity can be shown more clearly by the direction of Poynting vectors during the Dynamic process. From Fig. 7(a)–(d), we show the direction of Poynting vectors in the moments of Fig. 5(c)–(f), respectively. In Fig. 7, we see that there is an “intensity front” (shown by red dashed curve) which separates the two regions of the CS. At the right side of the front, the field intensity in the CS is much weaker than the outside and the Poynting vector directions are not regular (especially near the front). However, at the left-side region which is swept by the intensity front, the Poynting vectors are regular and nearly along the “cloaking rays” which was predicted at the coordinate transformation [40]. As the cloaking effect can be interpreted by the mimic picture that the light runs around the cloaking area through these curved cloaking rays, it is not surprising to find that the stable cloaking state is achieved when the intensity front sweeps through the whole CS and these optical rays are well constructed. The surprising thing is that the stable cloaking state seems to be constructed locally. We believe this property related to the original cloaking recipe [40], which makes the cloaking material almost impedance matched layer by layer. This also explains why there is no oscillation in the cloaking dynamic process generally. This picture can also interpret the strong scattered field in the dynamic process, as these “irregular rays” at the right-side region of the intensity front must be scattered strongly. Furthermore, we can use this picture to analyze the dynamic process of other incident waves, including the Gaussian beams, which are composed of different plane-wave components.

With this understanding, now we are ready to find the correlation between the relaxation time and the CS dispersion. It is well known that the field intensity (or energy) propagates at the group velocity V_g, which is controlled by the material dispersion. Therefore, the intensity front, which determines the dynamic process, should move in V_g. Thus, we can explain the results in Fig. 6, as our modification of the dispersive parameters can cause the V_g to change. However, because the cloaking material is a strong anisotropic material, the V_g at different directions could be very different. Can we predict more precisely which component dominates the relaxation time? The answer is “yes”. In Fig. 7(d), we can see that the stable energy flow in the CS is nearly along the θ direction at most regions of the CS. Then it is reasonable for us to argue that it is the component along the θ direction V_{g\theta }, not the component along the r direction V_{gr}, that dominates the relaxation time and the total scattered energy in the Dynamic process.

For the anisotropic cloaking material, the V_{g\theta } and V_{gr} can be expressed as: V_{g\theta }={[{\nabla }_{\mathit{k}}\omega (\mathit{k})]}_{\theta }=\left(\frac{2c}{\sqrt{{\epsilon }_z{\mu }_r}}\right)/\left(2+\frac{\omega }{{\epsilon }_z}\frac{\mathrm{d}{\epsilon }_z}{\mathrm{d}\omega }+\frac{\omega }{{\mu }_r}\frac{\mathrm{d}{\mu }_r}{\mathrm{d}\omega }\right)and V_{gr}={[{\nabla }_{\mathit{k}}\omega (\mathit{k})]}_r=\left(\frac{2c}{\sqrt{{\epsilon }_z{\mu }_{\theta }}}\right)/\left(2+\frac{\omega }{{\epsilon }_z}\frac{\mathrm{d}{\epsilon }_z}{\mathrm{d}\omega }+\frac{\omega }{{\mu }_{\theta }}\frac{\mathrm{d}{\mu }_{\theta }}{\mathrm{d}\omega }\right), where c is the velocity of light in vacuum.

In order to illustrate our prediction, the \sigma (t) under different V_{g\theta } and V_{gr} are investigated, respectively. First, we keep the V_{gr} unvaried by holding A_{\theta }, A_z and γ constant (keep \mathrm{d}{\epsilon }_z/\mathrm{d}\omega and \mathrm{d}{\mu }_{\theta }/\mathrm{d}\omega unchanged), and modify only A_r to change the V_{g\theta }. The results are shown in Fig. 6(c); when A_r is closer to one, the V_{g\theta } becomes smaller (with larger \mathrm{d}{\mu }_r/\mathrm{d}\omega), the relaxation time is longer and more energy is scattered in the dynamic process. Therefore, larger V_{g\theta } means better cloaking effect in the dynamic process. On the other hand, when we keep the V_{g\theta } unvaried and change V_{gr} by holding A_r, A_z, and γ constant and modifying A_{\theta }, the results are shown in Fig. 6(d). We find that the relaxation time is almost unchanged with the change of V_{gr}. Obviously, V_{g\theta } is the dominant element in the dynamic process. This conclusion can help us to design a better CS to defend the pulsive radars. In the expression of V_{g\theta }, it is also shown how to tune V_{g\theta } by modifying dispersion parameters.

It seems that the larger the V_{g\theta }, the better the cloaking effect will come out in the dynamic process. However, since V_g (and its components) cannot exceed c generally, there is a minimum limit for the relaxation time of the cloaking dynamic process. We can estimate it through dividing the mean length of the propagation rays by V_{g\theta }. In our model, the mean length is \text{π}(R_2+R_1)/2, about three wavelengths. So the relaxation time cannot be shorter than three periods. Figure 6 shows that our estimation is coincident with our simulation results. Actually, here we are facing a basic conflict to make a “better” CS, which is discussed more in our other works [9, 53]. The conflict lies in the fact that a pretty strong dispersion is required to realize a good stable cloaking effect at a certain frequency [9, 40], but in this research we show that the weaker dispersion can realize a better cloaking effect in the dynamic process. In the real design of the CS, there should be an optimized trade-off.

Based on causality, the limitation of the electromagnetic cloak with dispersive material is investigated in this section. The results show that perfect invisibility cannot be achieved because of the dilemma that either the group velocity V_g diverges or a strong absorption is imposed on the cloaking material. It is an intrinsic conflict which originates from the demand of causality. However, the total cross section can really be reduced through the approach of coordinate transformation. A simulation of finite-difference time-domain method is performed to validate our analysis.

4 Limitation of the electromagnetic cloak with dispersive material

Through the ages, people have dreamed of having a magic cloak whose owner cannot be seen by others. For this fantastic dream, plenty of work has been done by scientists all over the world. For example, researchers diminished the scattering or the reflection from objects by absorbing screens [54] and small, non-absorbing, compound ellipsoids [55]. More recently, based on the coordinate transformation, Pendry et al. theoretically proposed a general recipe for designing an electromagnetic cloak to hide an object from electromagnetic (EM) waves [40]. An arbitrary object may be hidden because it remains untouched by external radiation. Meanwhile, Ulf Leonhardt described a similar method where the Helmholtz equation is transformed to produce similar effects in the geometric limit [39, 47]. Soon, Cummer et al. simulated numerically(COMSOL) the cylindrical version of this cloak structure using ideal and nonideal (but physically realizable) electromagnetic parameters [88]. Especially, Schurig et al. experimentally demonstrated such a cloak by split-ring resonators [44]. In addition, Cai et al. proposed an electromagnetic cloak using high-order transformation to create smooth moduli at the outer interface and presented a design of a non-magnetic cloak operating at optical frequencies [42, 43]. According to the general recipe, the electromagnetic cloak is supposed to be perfect or “fully functioned” at a certain frequency as long as we can get very close to the ideal design although there is a singularity in the distribution, which has been elucidated further in several literatures [56, 95]. However, in all these pioneering works, the interests are mainly focused on single-frequency EM waves, so that the effects of the dispersion, which is related with basic physical laws, are not well studied. Suppose the dispersion is introduced into the study, can we have a deeper insight into the cloaking physics?

In this section, we will show that the ideal cloaking cannot be achieved due to another more basic physical limitation — the causality limitation (Basing on the same limitation, Chen et al. obtained a constraint of the bandwidth that limits the design of an invisibility cloak [57].). Starting from the dispersion relation and combining it with the demand of causality, we will demonstrate that the ideal cloaking will lead to the dilemma that either the group velocity V_g diverges or a strong absorption is imposed on the cloaking material. Our derivation and numerical experiments based on the FDTD methods will show that the absorption cross section will be pretty large and dominate the total cross section for a dispersive cloak, even with very small imaginary parts of permittivity and permeability.

Let us consider a more general coordinate transformation on an initial homogeneous medium with {\epsilon }_i={\mu }_i in r space: r^{\prime }=f(r), {\theta }^{\prime }=\theta, {\phi }^{\prime }=\phi. Following the approach in Refs. [58, 60], we get the following radius-dependent, anisotropic relative permittivity and permeability: {\epsilon }_{r^{\prime }}={\mu }_{r^{\prime }}={\epsilon }_i{\left[r/f(r)\right]}^2\mathrm{d}f(r)/\mathrm{d}r,{\epsilon }_{{\theta }^{\prime }}={\mu }_{{\theta }^{\prime }}={\epsilon }_i/[\mathrm{d}f(r)/\mathrm{d}r] and {\epsilon }_{{\phi }^{\prime }}={\mu }_{{\phi }^{\prime }}={\epsilon }_i/[\mathrm{d}f(r)/\mathrm{d}r]. We emphasize that since the transformation is directly acted on the Maxwell equations, the above equations are also suited for the imaginary parts of constitutive parameters, and all the physical properties of wave propagation in r space should be inherited in r^{\prime } space, such as the absorption. This is important for us to have consistent physical pictures in both spaces. At working frequency {\omega }_0, for a propagating mode with k-vector as \{k_{r^{\prime }},k_{{\theta }^{\prime }},k_{{\phi }^{\prime }}\} inside the cloak, we have the dispersion relation of the anisotropic material [58] as: k_{r^{\prime }}^2/n_{r^{\prime }}^2+k_{t^{\prime }}^2/n_{t^{\prime }}^2={\omega }^2/c^2, where k_{t^{\prime }}^2=k_{{\theta }^{\prime }}^2+k_{{\phi }^{\prime }}^2, n_{r^{\prime }}=\sqrt{{\epsilon }_{{\phi }^{\prime }}{\mu }_{{\theta }^{\prime }}}=n_i/[\mathrm{d}f(r)/\mathrm{d}r], n_{t^{\prime }}=n_{{\theta }^{\prime }}=n_{{\phi }^{\prime }}=\sqrt{{\epsilon }_{r^{\prime }}{\mu }_{{\theta }^{\prime }}}=n_{i}r/f(r), and n_i=\sqrt{{\epsilon }_i{\mu }_i}. Then we can define k_{r^{\prime }}=\omega /c{n}_{r^{\prime }}\cos \alpha and k_{t^{\prime }}=\omega /c{n}_{t^{\prime }}\sin \alpha, the group velocity can be obtained as:

\begin{array}{}\text{(4)}& V_g=c\sqrt{\frac{{(\cos \alpha )}^2}{n_{r^{\prime }}^2}+\frac{{(\sin \alpha )}^2}{n_{t^{\prime }}^2}}/[{(\cos \alpha )}^2\frac{m_{r^{\prime }}}{n_{r^{\prime }}}+{(\sin \alpha )}^2\frac{m_{t^{\prime }}}{n_{t^{\prime }}}]\end{array}

where m_{r^{\prime }}=n_{r^{\prime }}+\omega \mathrm{d}{n}_{r^{\prime }}/\mathrm{d}\omega and m_{t^{\prime }}=n_{t^{\prime }}+\omega \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega.

If the transformation has the following characteristic: f(r=0)=R_1, f(r=R_2)=R_2, then when r^{\prime }\to R_1 (or r\to 0), n_{t^{\prime }} will tend to zero, and the group velocity is approximated as:

\begin{array}{}\text{(5)}& V_g\approx \frac{c}{|\sin \alpha |\omega \frac{\mathrm{d}{n}_{t^{\prime }}}{\mathrm{d}\omega }}\end{array}

We will discuss Eq. (5) in two cases. The first case is with the finite \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega. Obviously, V_g will diverge when \sin \alpha \to 0 for any finite \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega. Such divergence is shown in Fig. 8 for a concrete example, in which the transformation is r^{\prime }=f(r)=(R_2-R_1)r/R_2+R_1 as Ref. [40], R_2=2R_1, thus n_{r^{\prime }}=2 and n_{t^{\prime }}=2-4/(r/R_1+2). The dispersion parameters are set as m_{r^{\prime }}=2.5, \omega \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega =4 at a working frequency. In Fig. 8, the curves of V_g vs. α are plotted for different R_1/r values. We can see that, for large R_1/r(r\to 0), the group velocity (more precisely, the tangential component of V_g) will diverge at both peaks around \alpha =0.

Because of the causality limitation, it is well-known that the group velocity cannot exceed c except in the “strong dispersion” frequency range (or called “resonant range”). However, if the working frequency is in the “strong-dispersion” range of the cloaking material, the absorption must be very strong and it will destroy the ideal cloaking obviously. Therefore, perfect invisibility cannot be achieved for the finite \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega because it will lead to superluminal velocity or strong absorption.

In addition, the curves with the criterion condition V_g=c on the plane [R_1/r, α] are plotted for different \omega \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega in Fig. 9. The region to the left of curves corresponds to and the region to the right corresponds to V_g>c. There exists a maximum \max \{R_1/r\} for each curve in order that V_g⩽c can hold for all the α. Especially, for the no-dispersion case \omega \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega =0, we can see that V_g>c at all R_1/r for large α values, which means the whole cloak is not physical if there is no dispersion. This “dispersion-is-required” conclusion can be generally derived from Eq. (4), and it is consistent with the analysis in Ref. [40]. From Fig. 9, we know that the larger the \omega \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega is, the larger the \max \{R_1/r\} will be. But anyway, for arbitrary finite \omega \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega, \max \{R_1/r\} cannot be infinite, so that the superluminal range always exists.

The second case of Eq. (5) is with divergent \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega. From the previous discussion, we know that if the ideal cloak exists, the cloak must be dispersive and \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega must be divergent when r\to 0. Actually, when r\to 0, since \sqrt{{\epsilon }_{r^{\prime }}}\propto r, \mathrm{d}{n}_{t^{\prime }}/\mathrm{d}\omega \propto (\mathrm{d}{\epsilon }_{r^{\prime }}/\mathrm{d}\omega )/\sqrt{{\epsilon }_{r^{\prime }}}is really divergent for non-zero \mathrm{d}{\epsilon }_{r^{\prime }}/\mathrm{d}\omega. From Eq. (5), we can see that now the V_g\to 0 for a finite \mathrm{d}{\epsilon }_{r^{\prime }}/\mathrm{d}\omega (generally true) at all α values except \alpha =0 (or π), so that the group velocity difficulty seems to be overcome. However, because of the causality limitation, the non-zero \mathrm{d}{\epsilon }_{r^{\prime }}/\mathrm{d}\omega means non-zero imaginary part of permittivity (non-zero dissipation). The non-zero dissipation and the almost-zero group velocity will result in very strong absorption. This means that the energy of rays near the inner cloaking radius R₁ is almost totally absorbed by the cloaking material. As we pointed out in the beginning the absorption in r^{\prime } space should also appear in r space, because of the consistency between two spaces, the strong absorption in r space can be interpreted in the following way. From the transformation (which is also suited for imaginary part), we can find that when r\to 0, the finite imaginary part in r^{\prime } space corresponds to the infinite imaginary part in r space, which also means very strong absorption in the initial homogeneous medium. Therefore, the perfect cloaking is still impossible because of the strong absorption which is enforced by the causality limitation.

For a 2-D coordinate transformation: r^{\prime }=f(r), {\theta }^{\prime }=\theta, z^{\prime }=z, the same conclusions of the causality limitation can be obtained through the similar analysis, although the coordinate transformation and the singularities are different from the three-dimensional case.

Next, we will discuss the physical meaning of the dilemma that either the group velocity V_g diverges or a strong absorption is imposed on the cloaking material. First, it is an intrinsic conflict that cannot be solved by any method, for example, “the system is imbedded in a medium” [40]. We believe that the ideal cloaking is impossible because of the causality limitation, and this conclusion is consistent with the statement of the previous studies [60] that the perfect invisibility is unachievable because of the wave nature of light. Second, we have to face the question “why is the causality violated for ideal cloaking which is based on the simple coordinate transformation?” Our answer is that the causality is guaranteed only by the Lorentz co-variant transformation, but the coordinate transformation for ideal cloaking is not Lorentz co-variant. Such violation is obvious if we suppose the initial medium in the r space is not dispersive, such as the vacuum, but as we have pointed out (also mentioned in Ref. [40]), the cloaking material (in r^{\prime } space) must be dispersive to avoid the group velocity over c. Such Lorentz co-variant violation is generally true for transformation optics as material parameters are non-relativistic, so the causality limitation should be checked widely. Third, from Eq. (4), we can find that not only the inner layers of the cloak (r^{\prime }\to R_1) but also the other layers (r^{\prime }>R_1) must be dispersive. For every layer, a certain dispersive strength is needed to avoid V_g>c.

In the following, we will validate that the total cross section can be reduced drastically, and that the perfect cloaking cannot be achieved because of strong absorption by FDTD numerical experiments. Compared with other frequency-domain simulation methods, such as the finite element methods or the transfer-matrix methods, the FDTD simulation can better reflect the real physical process of cloaking. For example, we note that the FDTD calculation will be numerically unstable when the dispersion is not included in the cloak’s material. For simplicity, the simulation is limited to 2-D cloak [88]. Without loss of generality, only TE modes are investigated in this study. (TE modes have the electric field perpendicular to the 2-D plane of our model.) Thus, the constitutive parameters involved here are {\epsilon }_{z^{\prime }},{\mu }_{r^{\prime }},{\mu }_{{\theta }^{\prime }}. The dispersion is introduced into our FDTD by standard Lorentz model:

\begin{array}{}\text{(6)}& \begin{array}{c}{\tilde{\epsilon }}_{z^{\prime }}(r^{\prime },\omega )=1+F_{z^{\prime }}(r^{\prime }){\omega }_{p{z}^{\prime }}^2/[{\omega }_{a{z}^{\prime }}{(r^{\prime })}^2-{\omega }^2-\mathrm{i}\omega {\gamma }_{z^{\prime }}]\hfill \\ {\tilde{\mu }}_{r^{\prime }}(r^{\prime },\omega )=1+F_{r^{\prime }}(r^{\prime }){\omega }_{p{r}^{\prime }}^2/[{\omega }_{a{r}^{\prime }}{(r^{\prime })}^2-{\omega }^2-\mathrm{i}\omega {\gamma }_{r^{\prime }}]\hfill \\ {\tilde{\mu }}_{{\theta }^{\prime }}(r^{\prime },\omega )=1+F_{{\theta }^{\prime }}(r^{\prime }){\omega }_{p{\theta }^{\prime }}^2/[{\omega }_{a{\theta }^{\prime }}{(r^{\prime })}^2-{\omega }^2-\mathrm{i}\omega {\gamma }_{{\theta }^{\prime }}]\hfill \end{array}\end{array}

\begin{array}{}\text{(7)}& r^{\prime }=\frac{R_2(\theta )-R_1(\theta )}{R_2(\theta )}r+R_1(\theta )\end{array}

\begin{array}{}\text{(8)}& {\theta }^{\prime }=\theta \end{array}

where {\omega }_{p{z}^{\prime }}, {\omega }_{p{r}^{\prime }}, {\omega }_{p{\theta }^{\prime }} are plasma frequencies, {\omega }_{a{z}^{\prime }}, {\omega }_{a{r}^{\prime }}, {\omega }_{a{\theta }^{\prime }} are atom resonated frequencies, {\gamma }_{z^{\prime }}, {\gamma }_{r^{\prime }}, {\gamma }_{{\theta }^{\prime }} are damping factors and F_{p{z}^{\prime }}, F_{p{r}^{\prime }}, F_{p{\theta }^{\prime }} are filling factors. In our simulation, an E-polarized time-harmonic uniform plane wave whose wavelength {\lambda }_0 in vacuum is 3.75 cm is incident from left to right. The real parts of the constitutive parameters at {\omega }_0=2\text{⁢π}c/{\lambda }_0 satisfy the cloaking coordinate transformation [61, 88], and they are {\mu }_{r^{\prime }}=(r^{\prime }-R_1)/r^{\prime }, {\mu }_{{\theta }^{\prime }}=1/{\mu }_{r^{\prime }}, {\epsilon }_{z^{\prime }}={[R_2/(R_2-R_1)]}^2(r^{\prime }-R_1)/r^{\prime }, where R_1is 0.665{\lambda }_0, R_2 is 1.33{\lambda }_0. And the dispersive parameters are set as follows: if {\epsilon }_{z^{\prime }}>1 then {\omega }_{a{z}^{\prime }}=1.4{\omega }_0, else {\omega }_{a{z}^{\prime }}=0.6{\omega }_0, {\omega }_{a{r}^{\prime }}=0.6{\omega }_0, {\omega }_{a{\theta }^{\prime }}=1.4{\omega }_0, {\gamma }_{z^{\prime }}={\gamma }_{r^{\prime }}={\gamma }_{{\theta }^{\prime }}={\omega }_0/100. {\omega }_{p{z}^{\prime }}={\omega }_{p{r}^{\prime }}={\omega }_{p{\theta }^{\prime }}=4{\omega }_0, F_{z^{\prime }}(r)=({\epsilon }_{z^{\prime }}-1)\frac{{({\omega }_{a{z}^{\prime }}^2-{\omega }_0^2)}^2+{\omega }_0^2{\gamma }_{z^{\prime }}^2}{({\omega }_{a{z}^{\prime }}^2-{\omega }_0^2){\omega }_{p{z}^{\prime }}^2}, F_{r^{\prime }}(r)=({\mu }_{r^{\prime }}-1)\frac{{({\omega }_{a{r}^{\prime }}^2-{\omega }_0^2)}^2+{\omega }_0^2{\gamma }_{r^{\prime }}^2}{({\omega }_{a{r}^{\prime }}^2-{\omega }_0^2){\omega }_{p{r}^{\prime }}^2}, F_{{\theta }^{\prime }}(r)=({\mu }_{{\theta }^{\prime }}-1)\frac{{({\omega }_{a{\theta }^{\prime }}^2-{\omega }_0^2)}^2+{\omega }_0^2{\gamma }_{{\theta }^{\prime }}^2}{({\omega }_{a{\theta }^{\prime }}^2-{\omega }_0^2){\omega }_{p{\theta }^{\prime }}^2}. In fact, these parameters have many possible choices. The different groups of parameters correspond to different dynamic processes which we will discuss in another section [62].

Figure 10 shows the snapshots of the electric-field distribution in two cases: the cloak with the perfect electric conductor (PEC) at radius R₁ [Fig. 10(a)], and the naked PEC with radius R₁ [Fig. 10(b)]. Obviously, the cloak is very effective. Quantitatively, we calculate the absorption cross section and the scattering cross section of the cloak at the stable state, and they are 0.67{\lambda }_0 and 0.24{\lambda }_0 respectively, while the scattering cross section of the naked PEC is 3.14{\lambda }_0. Therefore, with dispersive cloak, the total cross section is three times smaller, and the absorption cross section dominates as we predicted. To emphasize the huge absorption of the cloak, we use a common homogeneous isotropic media, with \epsilon =\mu =1.1 but all other parameters are the same as the cloak, to replace the cloaking material. Then we find the absorption cross section is only 0.089{\lambda }_0 which is about one order smaller. The reason for the strong absorption has been discussed before.

Now we have a full view of the cloaking recipe based on the coordinate transformation. First, the cloaking material must be dispersive, and the strong absorption cannot be avoided because of the causality limitation. Thus, it is not perfectly invisible. Second, the scattering cross section of the dispersive cloak could be small, so the scattered field is weak. Although the ideal invisibility is impossible, the cloaking recipe has still a main advantage. The “strong-absorption” and “weak-scattering” property means that the cloak almost cannot be observed except from the forward direction so that such a cloak can well defend the passive radars which detect the perturbation of the original field. It is well-known that at the Rayleigh scattering case, where the radius of the scatterer is much smaller than the wavelength, the absorption cross section could be larger than the scattering cross section because of the diffraction. The cloaking can be thought as a giant Rayleigh scattering case, where the light rays are forced to “diffract” around the cloaked area.

In conclusion, the properties of the dispersive cloak are investigated, and the limitation of causality is revealed. Our study shows that the superluminal velocity or the strong absorption cannot be overcome because of the intrinsic conflict between the coordinate transformation to obtain the cloaking and the causality limitation. In addition, we validate the results using a numerical simulation which is performed in FDTD algorithm with physical parameters. The numerical experiments show that the absorption cross section is dominant and the scattering cross section can be reduced significantly. The study gives us a full view of the cloaking recipe based on the coordinate transformation, and will have further profound influence on the related topics.

5 Hyper-interface as the bridge between radiative wave and evanescent wave

Many new phenomena are observed at the interfaces between metamaterial and common dielectric material, including the negative refraction which is found at the left-handed material (LHM) surface [63-65]. More interestingly, the evanescent wave (EW) could be amplified at the LHM interface so that the super-resolution can be achieved [66-68]. Besides the LHM, there is another class of anisotropic metamaterial, so-called “hyperbolic medium” (HM), in which one of the diagonal permittivity tensor components is negative and results in a hyperbolic dispersion. For convenience, we call the interface between a HM and a common dielectric material “hyper-interface” (HI). Some surprising electromagnetic properties of HI have been intensively studied recently [69-75]. For instance, HI can convert the EW into the radiative wave (RW) so that the sub-wavelength information can be observed at far-field, which is called “hyper-lens” effect [69-71]. Very recently it was found that when HI is perpendicular to one asymptote of HM dispersion, abnormal omnidirectional transmission occurs [76]. Although some theoretical and experimental works [72-75] have demonstrated that the EW can be really converted into RW by HI of the layered cylindrical HM, a full theory involving the “material dispersion” (will be explained later) has not been given so far. For metamaterial systems, if without physical dispersion, some abnormal optical properties cannot be clearly explained and the dynamic study of wave propagation cannot be carried out [78]. Even more seriously, the causality violation because of the superluminal group velocity (v_g>c) in HM is pointed out [10], which makes the observed hyperlens effect doubtful. To solve these problems and predict new phenomena, more robust theory with stronger base is needed. On the other hand, to compress and to stop (slow) light pulses are very essential to modern optical/photonics research and signal processing. Hence, a new mechanism, which can compress and stop (slow) light pulses and is frequency and direction insensitive, would induce wide interest in related directions.

In this section, we theoretically and numerically investigate the novel optical properties of flat HI [80-82], in which, unlike the cylindrical HI, the translational symmetry guarantees the simple physical picture for intuitive understanding, the quantitative study of the conversion between EW and RW etc. A general theory of HI is constructed with physical dispersion of HM. On the HI, not only the conversion from EW to RW (CER) of the hyperlens is confirmed, when RW is incident from HM to dielectric (the inverse process of hyperlens), but also the almost total conversion from RW to EW (CRE) can occur, i.e., there is “no-transmission and no-refelection” (NTNR). More importantly, we find that this is a new mechanism to compress and stop (slow) light pulses in wide frequency and direction range with many potential applications. Theoretically and numerically we demonstrate that the superluminal group velocity in hyperlens is artificial since the HM material dispersion is neglected in the previous study [10].

Finally, the feasibility to realize these functions on real structures is discussed. All theoretical results are demonstrated by FDTD simulations. Our model is as follows. Assuming two plane waves are incident to HI from HM and isotropic dielectric, and scattered from HI, as shown in the upper-right insert of Fig. 11. The HI is in the x–z plane, while the incident surface and both HM optical axes lie in the x–y plane. The incident waves are chosen as TM wave with field components (E_x,E_y,H_z) and same “parallel wave-vector” k_x. The HM is with the permittivity tensor as:

\begin{array}{}\text{(9)}& {\widehat{\epsilon }}_p=\left(\begin{array}{ccc}{\epsilon }_1(\omega )& 0& \\ 0& {\epsilon }_2& \end{array}\right)\end{array}

in its principle axes coordinate, where and {\epsilon }_2>0 in our assumption. Also, the permittivity of isotropic dielectric material is \epsilon. The essential point of our model is that the negative diagonal component is dispersive {\epsilon }_1={\epsilon }_1(\omega ), which is called material dispersion in our study. It is well known that dispersion is physically required for real metamaterials with abnormal effective constitutive coefficients, such as negative permittivity. We will see that the material dispersion will help us to obtain self-consistent explanation of abnormal optical properties of HI and to avoid causality violation.

For simplicity, the HM is nonmagnetic and Gaussian unit is employed throughout the paper. We define the angle between the HI (or x axis) and the positive-\epsilon principle axis of HM as θ, then the most general frequency contour in k space of HM is

\begin{array}{}\text{(10)}& \frac{{(k_x\cos \theta -k_y\sin \theta )}^2}{{\epsilon }_2}-\frac{{(k_y\cos \theta +k_x\sin \theta )}^2}{|{\epsilon }_1|}={\left(\frac{\omega }{c}\right)}^2\end{array}

In the HM region and the isotropic dielectric region the fields can be expressed uniformly as {\mathit{H}}^{\sigma }={\mathit{e}}_{\mathit{z}}(H_{iz}^{\sigma }+H_{sz}^{\sigma }){\mathrm{e}}^{\mathrm{i}{k}_{x}x-\mathrm{i}\omega t}, {\mathit{E}}^{\sigma }=[{\mathit{e}}_{\mathit{x}}(c_{ix}^{\sigma }{H}_{iz}^{\sigma }+c_{sx}^{\sigma }{H}_{sz}^{\sigma })+{\mathit{e}}_y(c_{iy}^{\sigma }{H}_{iz}^{\sigma }+c_{sy}^{\sigma }{H}_{sz}^{\sigma })]{\mathrm{e}}^{\mathrm{i}{k}_{x}x-\mathrm{i}\omega t} where the coefficients are defined as c_{\nu \xi }^{\sigma }=E_{\nu \xi }^{\sigma }/H_{\nu z}^{\sigma } with \nu =i,s for the incident fields to HI or the scattered fields from HI; \xi =x,y; \sigma =h,d for the HM region or dielectric region. Because of the translation symmetry of flat HI, the wave-vector parallel component k_x is continuous at both regions. In the following discussion, the k vectors are normalized by {\tilde{k}}_{\nu \xi }^{\sigma }=k_{\nu \xi }^{\sigma }/k_0 where k_0=\omega /c.

To explore the transmission and reflection properties of HI, we define a scattering matrix associated with the incident fields and out-going fields at HI,

\begin{array}{}\text{(11)}& {(H_{sz}^h,H_{sz}^d)}^{\mathrm{T}}=S_{2\times 2}{(H_{iz}^h,H_{iz}^d)}^{\mathrm{T}}\end{array}

where the superscript T means the matrix transpose, S_{11}=-(c_{ix}^h\epsilon +{\tilde{k}}_y^d)/(c_{sx}^h\epsilon +{\tilde{k}}_y^d), S_{12}=(-c_{ix}^h\epsilon +{\tilde{k}}_{iy}^d)/(c_{sx}^h\epsilon -{\tilde{k}}_{iy}^d), S_{21}=(c_{sx}^h-c_{ix}^h)\epsilon /(c_{sx}^h\epsilon +{\tilde{k}}_y^d), and S_{22}=-2{\tilde{k}}_{iy}^d/(c_{sx}^h\epsilon -{\tilde{k}}_{iy}^d). From the standard boundary conditions, the coefficients c_{i(s)x}^h,c_{i(s)x}^d are worked out to be c_{i(s)x}^h=(-{\tilde{k}}_{i(s)y}^h{\alpha }^2+{\tilde{k}}_x\alpha \gamma )/(|{\epsilon }_1|{\epsilon }_2);c_{i(s)x}^d=-{\tilde{k}}_{i(s)y}^d/\epsilon, where factors \alpha ,\gamma are defined as \alpha ={({\epsilon }_1{\sin }^2\theta -{\epsilon }_2{\cos }^2\theta )}^{\frac{1}{2}};\gamma =(|{\epsilon }_1|+{\epsilon }_2)\sin 2\theta /2\alpha, and the values of {\tilde{k}}_{i(s)y}^d=\pm \sqrt{\epsilon -{\tilde{k}}_x^2} and {\tilde{k}}_{i(s)y}^h and {\tilde{k}}_{i(s)y}^h=\{\pm {[|{\epsilon }_1|{\epsilon }_2({\tilde{k}}_x^2+{\alpha }^2)]}^{1/2}+{\tilde{k}}_x\alpha \gamma \}/{\alpha }^2 are uniquely determined by Eq. (10), respectively. From Eq. (11) we can easily get the reflection and transmission coefficients across the HI from up to down, or inverse. For the case of the hyperlens where the wave is incident from the isotropic medium to the HM, t_{dh}=S_{22}; r_{dh}=S_{12}. When {\tilde{k}}_x^2>\epsilon, the incident and the reflected waves in the isotropic dielectric are EWs with y-component wave-vectors as {\tilde{k}}_{iy}^d=\mathrm{i}\sqrt{{\tilde{k}}_x^2-\epsilon }=-{\tilde{k}}_{sy}^d\equiv \mathrm{i}\kappa. We note that, although single EW cannot carry net energy current (time averaged), two EWs, i.e., the incident and reflected EWs, can carry net energy current {\mathit{S}}_{iy} in isotropic dielectric medium, since the reflected EW gains an extra-phase from complex reflecting coefficient r_{dh}. The energy current {\mathit{S}}_{iy} carried by the two EWs can be converted by HI into the RW energy current {\mathit{S}}_{ty} in the HM:

\begin{array}{}\text{(12)}& \begin{array}{cc}|{\mathit{S}}_{iy}|\hfill & =\frac{\kappa }{\epsilon }\text{Im}{r}_{dh}=-\frac{2c_{sx}^h{\kappa }^2}{c_{sx}^{h\mathrm{\hspace{0.17em}2}}+{\kappa }^2}\hfill \\ \hfill & =|{\mathit{S}}_{ty}|=-\frac{1}{2}{|t|}^2{c}_{tx}\hfill \end{array}\end{array}

From Eq. (4), the hyper-lens effect and the image-improving by CER could be quantitatively studied.

After confirming CER on HI, it is natural to wonder if CRE can occur too, or if there are other novel phenomena on HI. Next we will study the inverse process of the hyperlens, i.e., the RW is incident from HM and the transmitted field is in the dielectric. For such inverse processes, there is a critical condition \theta ={\theta }_{\mathrm{c}}\equiv \arctan \sqrt{{\epsilon }_2/|{\epsilon }_1|}, which means HI (x axis) perpendicular to one of the hyperbola-dispersion asymptotes, or in other words, the asymptote is parallel with y axis now, as shown by the solid lines in Fig. 11. Under this critical condition, especially when the transmitted wave is EW, we will find CRE with NTNR, compressing and stopping light pulses, etc.

Before we get into detailed derivation, for the critical case (\theta ={\theta }_{\mathrm{c}}) we first present two seemingly conflicting conclusions of reflected wave from two different arguments, which will clearly show the most tricky point of HI.

The first argument is from the “intuitive way” which is based on Fig. 11. Since there is no reflection wave-vector on the dispersion curve to satisfy the k_x continuity, we intuitively expect that there should be no reflected wave with omnidirectional incidence. If the incident angle is as large as {\tilde{k}}_x^2>\epsilon then the transmitted field is EW, and since a single EW cannot carry the energy current, NTNR is the only possible choice and we expect that CRE will occur on HI. But from the second argument based on Eq. (4), we will obtain a different result. Since \theta ={\theta }_{\mathrm{c}} is a critical case, we should be more careful and discuss in a more subtle and strict way.

We first suppose \theta \ne {\theta }_{\mathrm{c}} as shown by dashed lines in Fig. 11, so the finite {\tilde{k}}_{ry} of reflected field for a fixed {\tilde{k}}_x can be found. Next we let the angle θ approach {\theta }_{\mathrm{c}} continuously (which can be realized physically by choosing a different direction of HI), then we find that {\tilde{k}}_{ry}\to \mathrm{\infty } when \theta \to {\theta }_{\text{c}} for a fixed {\tilde{k}}_x.

But surprisingly, when \theta \to {\theta }_{\text{c}}, the reflection coeffi-cient r_{dh}, calculated from Eq. (4) as {lim}_{\theta \to {\theta }_c}{r}_{dh}=[\epsilon {\tilde{k}}_x-\sqrt{|{\epsilon }_1|{\epsilon }_2(\epsilon -{\tilde{k}}_x^2)}]/[\epsilon {\tilde{k}}_x+\sqrt{|{\epsilon }_1|{\epsilon }_2(\epsilon -{\tilde{k}}_x^2)}] is not zero, and the reflected energy current is not zero either. Therefore, the theoretical result seems to be against our intuition.

To explain these conflicting results, we need to calculate the group velocity inside HM with material dispersion, which will also show that the superluminal group velocity is artificial. The most general expression of the group velocity of the reflected wave (which is also the group velocity of transmitted field in HM of the hyperlens case) can be obtained from Eq. (10) as:

\begin{array}{}\text{(13a)}& v_{gx}=\frac{2c{\epsilon }_1{\epsilon }_2(\alpha \gamma {\tilde{k}}_y^h+{\tilde{k}}_x{\beta }^2)}{{\epsilon }_2^2{\epsilon }_1{(\omega )}^{\prime }\omega {\tilde{k}}_y^{p2}-{\epsilon }_1^2{\epsilon }_2{(\omega )}^{\prime }\omega {\tilde{k}}_x^{p2}-2{\epsilon }_1^2{\epsilon }_2^2}\end{array}

\begin{array}{}\text{(13b)}& v_{gy}=\frac{2c{\epsilon }_1{\epsilon }_2(\alpha \gamma {\tilde{k}}_x-{\tilde{k}}_y^h{\alpha }^2)}{{\epsilon }_2^2{\epsilon }_1{(\omega )}^{\prime }\omega {\tilde{k}}_y^{p2}-{\epsilon }_1^2{\epsilon }_2{(\omega )}^{\prime }\omega {\tilde{k}}_x^{p2}-2{\epsilon }_1^2{\epsilon }_2^2}\end{array}

where {\tilde{k}}_x^p=({\tilde{k}}_x\cos \theta -{\tilde{k}}_y^h\sin \theta ) and {\tilde{k}}_y^p=({\tilde{k}}_y^h\cos \theta +{\tilde{k}}_x\sin \theta ) are the “k components” in the principal-axes coordinate of HM and \beta ={({\epsilon }_2{\sin }^2\theta -{\epsilon }_1{\cos }^2\theta )}^{1/2}. From Eq. (13), we find that, if the material dispersion of HM is neglected \partial {\epsilon }_1/\partial \omega =\partial {\epsilon }_2/\partial \omega =0, then we will obtain the superluminal group velocity as shown in Fig. 12(a). When \theta \to {\theta }_{\text{c}}, the v_g even diverges.

But with material dispersion, the x and y components of v_g is recalculated, and we find that there is no v_g>c in all cases, as shown in Fig. 12(b) in which two components of v_g versus \theta -{\theta }_{\mathrm{c}} based on Eq. (13), with the parameters {\epsilon }_2^{\prime }=0,{\epsilon }_2=1,{\epsilon }_1(\omega )=\left(1-{\omega }_p^2/{\omega }^2\right) [79] and \omega ={\omega }_p/\sqrt{2}, {\theta }_{\text{c}}=\text{π}/4. When approaching the critical angle \theta \to {\theta }_{\text{c}}, two components can be approximated as:

\begin{array}{}\text{(14a)}& \underset{\theta \to {\theta }_{\text{c}}}{lim}{v}_{gx}\sim \frac{1}{{\epsilon }_1^{\prime }(\omega )}{(k_{ry}^h)}^{-1}\end{array}

\begin{array}{}\text{(14b)}& \underset{\theta \to {\theta }_{\text{c}}}{lim}{v}_{gy}\sim \frac{k_x}{{\epsilon }_1^{\prime }(\omega )}{(k_{ry}^h)}^{-2}\end{array}

Since {\tilde{k}}_{ry}^h\to \mathrm{\infty } at the critical angle {\theta }_{\mathrm{c}}, the group velocity of the reflected wave should be zero {\mathit{v}}_g=0 at the critical angle, as shown in Fig. 12(b) too.

What does the zero-group velocity of reflected wave mean? The analysis will give us a clear answer. As we have pointed out, since the reflected energy current {\mathit{S}}_r is not zero and {\mathit{S}}_r=v_{g}W, where W is the energy density of reflected wave, the energy density W must be infinitely large at the critical angle. From Eq. (10), we can obtain that the electric field of reflected wave |{\mathit{E}}_{\mathit{r}}^{\mathit{h}}| is really divergent at the critical angle. The divergent field strength means that it needs infinite long time to accumulate energy at HI for the reflected field. In other words, there is no reflected wave physically, as our intuition has told us. When the incident angle is as large as {\tilde{k}}_x^2>\epsilon, since the energy of incident RW cannot be transmitted, and also cannot be reflected, the only answer is that the energy is stored at the HI or CRE occurs. Thus, we can have a self-consistent explanation for those seemingly conflicting results.

To confirm our theoretical discussion at {\theta }_{\mathrm{c}}, the FDTD simulation [51] with strict physical HM dispersion (Drude mode) which satisfies Kramar–Kronig relation, is done. The parameter of HM and dielectric are {\epsilon }_1=-3,{\epsilon }_2=3 , and {\epsilon }_l=1.1. For the case ({\tilde{k}}_x>\epsilon), as shown in Fig. 13(a), a light beam is incident from HM to HI in {45}^{\circ } angle. As we predicted, there is no reflection and no transmission, and the field energy is accumulated at HI and stopped there. More detailed observation shows that at the boundary the field energy is mainly at the dielectric side, as shown in Fig. 13(b).

We also checked the group velocities of hyperlens cases and find no violation of the causality. Actually, in FDTD simulation, if there is superluminal group velocity the program will be numerically unstable. The dynamic study, such as with the pulse incidence, can reveal some more interesting phenomena of HI. Since the group velocity along HI is zero at the NTNR case as discussed, we expect that the pulse energy will accumulate on HI and stay at the incident position until it is dissipated because of absorption of HM.

The numerical experiments with incident pulses by FDTD are also done. As shown in Fig. 13(c) and (d), two pulses arrive at the HI at different times, then they stop at the incident positions on HI. The pulse vertical length is compressed to almost zero, but their width stays the same so they are still well separated in Fig. 13(d). We emphasize here that this is a novel mechanism to compress and stop (slow) light pulses with special advantages. The first advantage is that this mechanism works at wide frequency and wide incident-angle range, which is confirmed by FDTD simulation in Fig. 13 with incident of pretty short pulses. The frequency and incident-angle insensitivity results from the fact that the mechanism possesses from a simple geometry property, i.e., the HI (x axis) perpendicular to one of the hyperbola-dispersion asymptotes. The second is that the decay (because of dissipation) of trapped field on HI is much slower than in common metallic material since the trapped field energy is mainly in the dielectric side as shown in Fig. 13(b). The third is that the trapped signals are easy to take out (read) as they are on the interface. Because of these advantages, HI could be used as a removable recorder (Dynamic memory) in optical/photonic signal processing, or as a wide-frequency wide-angle light trapper in photovoltaic devices.

However, we should point out that the above theoretical and numerical studies lie in the assumption of the ideal hyper-dispersion, which is still void when k_{ry}\to \mathrm{\infty }. In reality, such HM does not exist, therefore we need to study the limit of hyper-dispersion of realizable HM. HM can be realized by many structures, i.e., 1-D periodic metal-dielectric binary layers [83, 84] or 2-D periodic metallic lines [85], as shown in Fig. 14(a). For these structures, the dispersion relation can be calculated exactly. In Fig. 14, the calculated frequency contour of a 1-D metal-dielectric binary layers is shown, from which we can see that the effective HM medium is not available anymore when |k| approaches \text{π}/a. Based on this limit, we can roughly estimate the slow limit of group velocity by v_{gx}\sim 1/{\epsilon }_1^{\prime }(\omega ){(k_{ry}^h)}^{-1}\propto 1/{\gamma }_s and v_{gy}\sim k_x/{\epsilon }_1^{\prime }(\omega ){(k_{ry}^h)}^{-2}\propto 1/\gamma _s^{}{}_{}{}^2, where {\gamma }_s=k_{ry}/k_0 is the slowing coefficient. For the 2-D metallic-line structure, from the modern technical limit we assume the smallest lattice constant as a=10 nm. If the incident is the micro-wave \omega =5.8 GHz ({\gamma }_s\sim {10}^7) and {\epsilon }_1^{\prime }(\omega )=6.9\times {10}^{-10}s as in Ref. [86], we obtain v_{gx}\sim 4.6 m/s which means considerably slow light though not totally stopped, and v_{gy}\sim 7.07\times {10}^{-8}m/s which means that the strongly-compressed light pulses can be easily achieved.

In conclusion, we have theoretically and numerically investigated the optical properties of HI. The theory with dispersion of metamaterial is constructed and the hyperlens effect of CER is confirmed. At the inverse process of hyperlens, the abnormal phenomena of CRE with NTNR and a novel mechanism to compress and stop light in wide frequency range are revealed. Basing on the calculated group velocity, we demonstrate that the previously-pointed-out superluminal group velocity in HM is artificial since the material dispersion is neglected. FDTD simulations confirm that the HI has the potential to be a removable optical/photonic recorder, or a wide-frequency wide-angle light trapper. Finally, the realizability of these phenomena on the real metallic structures is discussed. Obviously, the new mechanism works not only for electromagnetic waves, but also for acoustic or matter waves if hyperbolic dispersion is available, which means that more interesting phenomena and applications are waiting for further theoretical and experimental research.

6 Remote control of light behavior by transformation optical devices

The dream of light-behavior remote control has haunted the human beings for thousands of years and there are legends with such imagination. In reality, although different methods of light-behavior controlling have been successfully used, the technique of light-behavior remote control has not come into being yet. One basic method to control light behavior in practice is using different cavities which are widely applied in electromagnetic/optical devices, such as resonators, modulators, switches and lasers [87]. Recently, a totally new direction of modern photonics/optics, i.e., the “transformation optics”, was setup by the research of an invisibility cloak [39, 40], which is designed from geometry transformation [39, 40, 42-44, 60, 88-98]. After that, the devices like super scatter [99] and anti-cloak [100] etc., were designed based on transformation optics to change light propagation. Furthermore, a new concept of “complementary medium” [102-104] is generated in transformation optics, i.e., a designed medium which can cancel the scattering effect of another existing medium in distance so that the scattering effect or the reflecting effect of “the complementary medium and the existing medium” is zero. However, the light-behavior remote control by transformation optics has not yet been put forward and further studied. In this paper, combining the cavity tuning methods and the theory of transformation optics, we propose a general light-behavior remote control method which uses a transformation optical device to remotely modify the cavity. The essential idea is that the transformation optical device is designed as the complementary medium of a part of optical cavity, e.g., the cavity mirror, then its reflection can be remotely tuned (even to zero-reflection), and the optical property of a cavity could be totally changed as expected and the light-behavior remote control can be realized. In this study, the remote cloak [101] is widely used to tune the cavity, but as we note, besides the remote cloak, other kinds of devices designed from more general folded geometry transformation [105] could be utilized for light-behavior remote control as well.

The method of light-behavior remote control in this work has some advantages. First, since the cavity is not modulated directly, there is no physical change of the cavity and no damage caused. Therefore, the cavity can be conveniently used repeatedly. Second, the modulating range can be larger. For example, it can modulate the quality factor Q of a cavity in several orders (even to almost zero) or tune the resonant frequency in a large range. Third, since the cavity is modulated at a distance and the quality factor Q could be tuned to a very low value, even to almost zero, there seems to be no cavity at all, and both the controller and the cavity could be well hidden.

To demonstrate that our method works, three basic schemes by tuning the cavity at a distance with a remote cloak are shown. The numerical simulations are in FDTD method and FEM. The novel phenomena, like “remotely enhancing the brightness of a cavity”, “hiding a laser which can be switched on remotely at anytime”, and “changing resonant-modes of a cavity in a large frequency range”, are shown to be theoretically possible.

6.1 The physical design and basic schemes

In this section, we will explain the basic design and three schemes of light-behavior remote control. Generally, for a cavity, the essential parameters are the quality factor Q and the resonant frequency {\omega }_0. We will show that both Q and {\omega }_0 could be remotely tuned. Our discussions are based on the Fabry–Perot (FP) cavity, since it is the simplest and most typical cavity that is widely used in optical/photonic devices. As the FP cavity shown in Fig. 15(a), there are two mirrors S₁ and S₂ at both sides of the  cavity  whose  reflection  coefficients are  r₁  and   r₂, respectively. Between the two mirrors there is a kind of filling medium in which the wave vector is k, thereby the cavity quality factor Q is determined by the value of r_1{r}_2 and the resonance frequency {\omega }_0 is determined by kd=n where nπ is an integer number.

If we can tune Q and {\omega }_0 of the FP cavity at a distance by changing the key parameters r_i (i=1,2) and kd with the remote cloak, the light-behavior remote control could be realized. Although this design is a general one, our model is a 2-D one and our study focuses on the E-polarization modes in which only the {\mu }_r, {\mu }_{\theta }, and {\epsilon }_z are considered, as shown in Fig. 15. For simplicity, r₁ is kept unchanged in our study, and only r₂ and kd are tuned. For the first two schemes in [Fig. 15(a) and (b)], the S₃ in the remote cloaks are designed as the complementary medium of the right mirror S₂, then the mirror can be totally canceled for optical effects. In other words, the reflection coefficient r₂ can be tuned to zero, and then the quality factor Q can be tuned considerably. Actually, with different remote cloak designs, Q can be tuned from the original value to zero, since S₃ can be designed to cancel only part of mirror S₂, then “physically effective thickness” of the mirror can be changed from the original thickness to zero. The difference between Fig. 15(a) and (b) is that Fig. 15(a) is with absorptive medium in the FP cavity, while Fig. 15(b) is with gain medium. In scheme 3 [Fig. 15(c)], the FP cavity which is made of a dielectric slab and S₃ is designed to cancel part of the slab, then the effective optical thickness kd of the cavity can be tuned and the resonant frequency {\omega }_0 is changed. The numerical simulation and detailed study will be shown in the next section.

The remote cloaks in the three schemes are similar to the structure in Ref. [101]. A remote cloak is composed of two parts which are shown in Fig. 15(a). Part I is the core area and part II is the cloaking shell, where a and b are the inner and the outer radii of the cloaking shell, respectively, and c is the effective radius of the remote cloak. In scheme 1 and scheme 2, S₃ inside the cloaking shell is designed as the complementary medium of mirror S₂, so r₂ of S₂ could be changed. In scheme 3, S₃ is designed to counteract part D_{\Delta } of the dielectric layer D_d, then the resonant frequency of D_d will be shifted. In order for these effects, we choose the permeability and the permittivity parameters as: {\mu }_{r,I}=1, {\mu }_{\theta ,I}=1, and {\epsilon }_{z,I}={\left(c/a\right)}^2 in the core area, {\mu }_{r,II}=\frac{r-a}{r}+\frac{c}{r}\frac{b-a}{b-c}, {\mu }_{\theta ,II}=1/{\mu }_{r,II}, and {\epsilon }_{z,II}={\left(\frac{b-c}{b-a}\right)}^2{\mu }_{r,II} in the cloaking shell region, and {\mu }_{r,S_3}={\mu }_{r,II}{\mu }_{r,i}, {\mu }_{\theta ,S_3}={\mu }_{\theta ,II}{\mu }_{\theta ,i}, and {\epsilon }_{z,S_3}={\epsilon }_{z,II}{\epsilon }_{z,i} for S₃, where {\mu }_{r,i}, {\mu }_{\theta ,i} and {\epsilon }_{z,i} are the permeability and the permittivity of the part of FP cavity to be counteracted.

6.2 Numerical simulations and discussion

To demonstrate the effect of our designs of light-behavior remote control, we have done numerical simulations of the three schemes.

In scheme 1, we illustrate how to remotely tune the output energy current (the brightness) from an absorptive FP cavity. For the FP cavity in Fig. 16(a), S₁ is a highly-reflective metal mirror (r_1\sim 1), while S₂ is a dielectric mirror ({\epsilon }_{z,s_2}>0) whose effect can be modified by a remote cloak (with the same dispersive model in Ref. [94]). Between the two mirrors, there is the absorptive medium, whose permittivity is \epsilon ={\epsilon }^{\prime }+\mathrm{i}{\epsilon }^{\prime \prime }, where {\epsilon }^{\prime }=1.01 and {\epsilon }^{\prime \prime }=0.1. A linear light source is set at the middle position of the cavity and the light emits from S₂. Also, we detect the output energy current on the boundary B₁. Without the remote cloak, the field by FDTD simulation is also shown in Fig. 16(a), which demonstrates that the light inside the FP cavity is reflected by the two mirrors and forms a resonant mode. Only a small part of light energy can emit from the mirror S₂ because of the absorption in the cavity. Next, we set a remote cloak on the right side to the FP cavity to modify the reflection coefficient r₂ of mirror S₂, as shown in Fig. 16(b). Since S₃ is the complementary medium of the mirror S₂, the effective thickness of S₂ can be reduced. When S₃ is designed to exactly cancel the mirror S₂ (the effective thickness of S₂ is zero), the light from the source inside the cavity propagates rightward like in a free space, as depicted by the field shown in Fig. 16(b). Since the light is not reflected between the two mirrors any more, much less energy is absorbed, so the output energy current can be larger than that of the original FP cavity.

We have also carried out a quantitative study of this effect. In Fig. 16(c), the output energy currents vs. imaginary part of dielectric constant {\epsilon }^{\prime \prime } are compared in two cases, case 1 with only an FP cavity, case 2 with an FP cavity and a remote cloak. Since the output energy current J is proportional to the decay factor \exp (-\alpha ) where \alpha \propto Q{\epsilon }^{\prime \prime }, the brightness of a cavity will exponentially depend on the quality factor Q and the material dissipation {\epsilon }^{\prime \prime }. The numerical results of \ln J vs. {\epsilon }^{\prime \prime } in Fig. 16 clearly show such dependence, and the \ln J of the case 2 with smaller Q since mirror S₂ is cloaked is obviously larger. (\ln J still decreases since the electromagnetic wave travels through the absorptive medium.)

In scheme 2, we demonstrate how to remotely control the lasing behaviors by the remote cloak to modify the mirror S₂. The setup of the FP cavity in Fig. 17 is the same as scheme 1, but with gain material inside, which is simulated by four-level atomic model [106]. To simulate the spontaneous emission, some quasi-monochromatic point sources are set randomly inside the cavity. According to the laser theory, when the gain is less than the threshold, there is only amplified spontaneous emission, but if the gain is over the threshold, the lasing will appear, which is shown in Fig. 17(a). The lasing threshold is proportional to the inverse of quality factor 1/Q. In Fig. 17(b), the remote cloak which can exactly cancel the mirror S₂ is added to the right side of the laser cavity. Since the light is not reflected by S₂ anymore and Q is very small, the threshold becomes so high that the lasing mode is not excited, and only the amplified spontaneous emitting field is observed. To display the difference, we also show the output energy current versus time for the two cases in Fig. 17(c). The blue dashed line depicts the lasing process in which J increases evidently (stable at about t/T=2500), while the red solid line shows a constant low output energy for the amplified spontaneous emission. Scheme 2 is not only a way to remotely switch off or on a laser, but also a good way to hide a prepared laser. If we cancel both the mirrors S₁ and S₂ of our model by two complementary mirrors in cloaks, then observers far away will not find any difference and such concealing could be a plot for future science fictions.

Besides the remote tuning of a mirror of FP cavity in scheme 1 and scheme 2, we can also tune the optical path kd to shift the resonant frequency, which is shown in scheme 3. As depicted in Fig. 15(c), the FP cavity in scheme 3 is a little different, only a thick dielectric layer, for which the two mirrors are the two interfaces between the dielectric and vacuum. If a remote cloak is added on the right side of a dielectric, then the “effective optical path” of the cavity becomes k{d}^{\prime }, where d^{\prime }=d-\Delta. Hence, the resonant condition becomes k{d}^{\prime }=n\text{π}. To demonstrate this effect, we calculate transmission spectrum in two cases by FEM, as shown in Fig. 18(a). In case 1, the spectrum of the original dielectric layer D_d is shown by the blue line marked with stars. We can see that there are three original resonant modes. In case 2, the spectrum of the FP cavity with a remote cloak on the right side is shown by the red dashed line marked with circles. From the spectrum, we can observe that the resonant modes shift. Numerically, we also have calculated the spectrum of a dielectric layer with thickness d-\Delta, shown by a black line in Fig. 18(a). We can see that these resonant modes of a shorter cavity are exactly the same with case 2, just as we expected. In Fig. 18(a), the largest shift range is about 15\%. In fact, with this method the resonant frequency can be tuned from a small shift to a very large one with different cloak design.

Furthermore, we can extend light-behavior remote control method to remotely tune some complex cavities. For example, we can tune 1-D photonic crystals (PhCs) by a transformation-optics device which is similar to that in Ref. [107]. The scheme is depicted in Fig. 18(b). By changing the effective optical path of every layer by transformation optical device, the photonic bandgap structure can be tuned in a wide range without any physical change PhCs. Similarly, the defect mode inside PhCs can also be tuned or generated by the transformation optical devices.

6.3 Conclusion

In conclusion, a light-behavior remote control method based on the transformation optics is proposed in this work. The tune of the quality factor Q and the resonant frequency {\omega }_0 of the cavities in wide ranges are demonstrated. By modifying these important characters of cavities, the light behavior can be remotely controlled without any physical change or damage to the cavities. In this work, we present three schemes, i.e., the output energy current of an absorptive cavity, the lasing behavior of a cavity with gain, and the cavity resonant frequency or the photonic band-gap of PhCs could be controlled by transformation optical devices. This work has proposed a new way for the application of transformation optical devices. With the fabrication breakthrough of metamaterial and transformation optical devices, we believe that the light-behavior remote control will be widely used, not only in optical/photonic devices, but also in electromagnetic devices, such as radar or antenna, etc.

7 Summary

In summary, in this paper we have investigated the Dynamic properties of the metamaterial with dispersion taken into account. Deeper physical pictures involved with metamaterial systems have been revealed, such as the temporal-coherence gain of the image formed by negative-index superlens with the quasi-monochromatic random source, the dynamic process of the dispersive cloak before achieving the stable state and the essential element to determine relaxation time and scattered energy, the limitation to the perfect invisibility of the electromagnetic cloak due to the material’s dispersion and the physical causality limit, the method of expanding the cloaking frequency range, etc. Furthermore, we also give two applications of the dynamic properties of the metamaterial system. One of them is stopping light on the interface between metamaterial and common medium, which is, to a large degree, dependent on the dispersion and anisotropy of the metamaterial. Another application is remotely controlling the quality factor Q and the resonant frequency {\omega }_0 of the cavity by a second transformation optical device, which is also based on the realization of metamaterial with special dependent anisotropy. Actually, we believe that only with the well-constructed dynamic picture is, a deeper understanding of the abnormal optical/photonic properties of metamaterials possible. All these research work show that the dynamic study of metamaterials can lead us to many new interesting topics, which are yet waiting for further research.

References

[1] V. G. Veselago, Sov. Phys. Usp., 1968, 10: 509

[2] D. R. Smith, W. J. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, Phys. Rev. Lett., 2000, 84: 4184

[3] R. A. Shelby, D. R. Smith, S. C. Nemat-Nasser, and S. Schultz, Appl. Phys. Lett., 2001, 78: 489

[4] R. A. Shelby, D. R. Smith, and S. Schultz, Science, 2001, 292: 77

[5] E. Cubukcu, K. Aydin, E. Ozbay, S. Foteinopoulou, and C. M. Soukoulis, Phys. Rev. Lett., 2003, 91: 207401

[6] E. Cubukcu, K. Aydin, E. Ozbay, S. Foteinopoulou, and C. M. Soukoulis, Nature (London), 2003, 423: 604

[7] P. Markos and C. M. Soukoulis, Phys. Rev. B, 2001, 65: 033401

[8] D. R. Smith and N. Kroll, Phys. Rev. Lett., 2000, 85: 2933

[9] P. Yao, Z. Liang, and X. Jiang, Appl. Phys. Lett., 2008, 92: 031111

[10] H. L. Luo, W. Hu, and Z. Z. Ren, Europhys. Lett., 2006, 74: 1081

[11] J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, IEEE Trans. Microwave Theory Tech., 1999, 47: 2075

[12] J. B. Pendry, Phys. Rev. Lett., 2000, 85: 3966

[13] J. B. Pendry, Phys. Rev. Lett., 2003, 91: 099701

[14] D. R. Smith, D. Schurig, M. Rosenbluth, S. Schultz, S. A. Ramakrishna, and J. B. Pendry, Appl. Phys. Lett., 2003, 82: 1506

[15] G. Gomez-Santos, Phys. Rev. Lett., 2003, 90: 077401

[16] S. Foteinopoulou, E. N. Economou, and C. M. Soukoulis, Phys. Rev. Lett., 2003, 90: 107402

[17] J. B. Pendry and D. R. Smith, Phys. Rev. Lett., 2003, 90: 029703

[18] X. S. Rao and C. K. Ong, Phys. Rev. B, 2003, 68: 113103

[19] X. S. Rao and C. K. Ong, Phys. Rev. E, 2003, 68: 067601

[20] M. W. Feise and Y. S. Kivshar, Phys. Lett. A, 2005, 334: 326

[21] L. Zhou and C. T. Chan, Appl. Phys. Lett., 2005, 86: 101104

[22] Y. Zhang, T. M. Grzegorczyk, and J. A. Kong, PIER, 2002, 35: 271

[23] L. Chen, S. He, and L. Shen, Phys. Rev. Lett., 2004, 92: 107404

[24] R. W. Ziolkowski and E. Heyman, Phys. Rev. E, 2001, 64: 056625

[25] S. A. Cummer, Appl. Phys. Lett., 2003, 82: 1503

[26] P. F. Loschialpo, D. L. Smith, D. W. Forester, F. J. Rachford, and J. Schelleng, Phys. Rev. E, 2003, 67: 025602(R)

[27] R. Merlin, Appl. Phys. Lett., 2004, 84: 1290

[28] C. Luo, S. G. Johnson, J. D. Joannopoulos, and J. B. Pendry, Phys. Rev. B, 2003, 68: 045115

[29] N. Engheta, IEEE Antennas and Wireless Propagation Lett., 2002, 1: 10

[30] I. V. Shadrivov, A. A. Sukhorukov, and Y. S. Kivshar, Phys. Rev. E, 2003, 67: 057602

[31] A. C. Peacock and N. G. R. Broderick, 2003, 11: 2502

[32] B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics, New York: John Wiley & Sons, 1991

[33] L. Mandel and E. Wolf, Optical Coherence and Quantum Optics, Cambridge: Cambridge University Press, 1995

[34] M. O. Scully and M. S. Zubairy, Quantum Optics, Cambridge: Cambridge University Press, 1997

[35] X. Jiang and C. M. Soukoulis, Phys. Rev. Lett., 2000, 85: 70

[36] In our source frequency range, the index range is about -1−0.0029i±(0.006+10⁻⁶i), so the focal length defference and reflection are very small.

[37] The “group velocity” is not a well-defined value if the working frequency ω₀ is near the resonant frequency ω_a of the NIM. But the GRT is still well-defined.

[38] X. Y. Jiang, et al., unpublished

[39] U. Leonhardt, Science, 2006, 312: 1777

[40] J. B. Pendry, D. Schurig, and D. R. Smith, Science, 2006, 312: 1780

[41] S. A. Cummer, B. I. Popa, D. Schurig, and D. R. Smith, Phys. Rev. E, 2006, 74: 036621

[42] W. Cai, U. K. Chettiar, A. V. Kildishev, and V. M. Shalaev, Nature Photonics, 2007, 1: 224

[43] W. Cai, U. K. Chettiar, A. V. Kildishev, V. M. Shalaev, and G. W. Milton, Appl. Phys. Lett., 2007, 91: 111105

[44] D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, Science, 2006, 314: 977

[45] H. Chen, B. I. Wu, B. Zhang, and J. A. Kong, Phys. Rev. Lett., 2007, 99: 063903

[46] Z. Ruan, M. Yan, C. W. Neff, and M. Qiu, Phys. Rev. Lett., 2007, 99: 113903

[47] U. Leonhardt, New J. Phys., 2006, 8: 118

[48] S. A. Cummer, Appl. Phys. Lett., 2003, 82: 2008

[49] X. Jiang, W. Han, P. Yao, and W. Li, Appl. Phys. Lett., 2006, 89: 221102

[50] P. Yao, W. Li, S. Feng, and X. Jiang, Opt. Express, 2006, 14: 12295

[51] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 2nd Ed., Boston: Artech House, 2000

[52] The approximation of permittivity ε and permeability μ in this section are: Re⁢[μθ⁢(r,ω0)]max=20, Re⁢[μr⁢(r,ω0)]min=1/20 and Re⁢[εz⁢(r,ω0)]min=1/5.

[53] H. Chen, Z. Liang, P. Yao, X. Jiang, H. Ma, and C. T. Chan, Phys. Rev. B, 2007, 76: 241104

[54] R. L. Fante and M. T. McCormack, IEEE Trans. Antennas Propag., 1968, 30: 1443

[55] M. Kerker, J. Opt. Soc. Am., 1975, 65: 376

[56] Z. C. Ruan, M. Yan, C. W. Neff, and M. Qiu, Phys. Rev. Lett., 2007, 99: 113903

[57] H. Y. Chen, X. Y. Jiang, and C. T. Chan, arXiv: 0707.1126v2, 2007

[58] D. Schurig, J. B. Pendry, and D. R. Smith, Opt. Express, 2006, 14: 9794

[59] H. Chen and C. T. Chan, Appl. Phys. Lett., 2007, 90: 241105

[60] U. Leonhardt and T. G. Philbin, New J. Phys., 2006, 8: 247

[61] When r′→R1, μθ′ will tend to infinite. In order that it can be realizable in our numerical simulation, we limit its maximum value to 103.

[62] Z. X. Liang, P. J. Yao, X. Y. Jiang, and X. W. Sun, unpublished

[63] R. A. Shelby, D. R. Smith, and S. Schultz, Science, 2001, 92: 792

[64] D. R. Smith and D. Schurig, Phys. Rev. Lett., 2003, 90: 077405

[65] C. G. Parazzoli, R. B. Greegor, K. Li, B. E. C. Koltenbah, and M. Tanielian, Phys. Rev. Lett., 2003, 90: 107401

[66] J. B. Pendry, Phys. Rev. Lett., 2000, 85: 3966

[67] N. Fang, H. Lee, C. Sun, and X. Zhang, Science, 2005, 308: 534

[68] T. Taubner, D. Korobkin, Y. Urzhumov, G. Shvets, and R. Hillenbrand, Science, 2006, 313: 1595

[69] J. B. Pendry, Opt. Express, 2003, 11: 755

[70] Z. W. Liu, H. Lee, Y. Xiong, C. Sun, and X. Zhang, Science, 2007, 315: 1686

[71] I. I. Smolyaninov, Y. J. Hung, and C. C. Davis, Science, 2007, 315: 1699

[72] Z. Jacob, L. V. Alekseyev, and E. Narimanov, Opt. Express, 2006, 14: 8247

[73] A. Salandrino and N. Engheta, Phys. Rev. B, 2006, 74: 075103

[74] Z. B. Jacob, L. V. Alekseyev, and E. Narimanov, J. Opt. Soc. Am. A, 2007, 24: A52

[75] J. B. Pendry and S. A. Ramakrishna, J. Phys.: Condens. Matter, 2002, 14: 8463

[76] X. Li, Z. X. Liang, X. H. Liu, X. Y. Jiang, and J. Zi, Appl. Phys. Lett., 2008, 93: 171111

[77] P. W. Milonni, Fast Light, Slow Light and Left-Handed Light, Institute of Physics, Great Britain: CRC Press, 2004

[78] As for the explanation of the static solution of the HI, it will be given in another paper.

[79] D. R. Smith and N. Kroll, Phys. Rev. Lett., 2000, 85: 2933

[80] V. A. Podolskiy and E. E. Narimanov, Phys. Rev. B, 2005, 71: 201101

[81] R. Wangberg, J. Elser, E. E. Narimanov, and V. A. Podolskiy, J. Opt. Soc. Am. B, 2006, 23: 498

[82] A. A. Govyadinov and V. A. Podolskiy, Phys. Rev. B, 2006, 73: 115108

[83] A. Salandrino and N. Engheta, Phys. Rev. B, 2006, 74: 115108

[84] N. Fang, H. Lee, C. Sun, and X. Zhang, Science, 2005, 308: 534

[85] J. Yao, Z. W. Liu, Y. M. Liu, Y. Wang, C. Sun, G. Bartal, A. M. Stacy, and X. Zhang, Science, 2008, 321: 930

[86] J. B. Pendry and A. J. Holden, and W. J. Stewart, Phys. Rev. Lett., 1996, 76: 4773

[87] B. E. A. Saleh and M. C. Teich, Foundamentals of Photonics, New York: John Wiley & Sons, Inc., 1991

[88] S. A. Cummer, B. I. Popa, D. Schurig, and D. R. Smith, Phys. Rev. E, 2006, 74: 036621

[89] S. Zhang, D. A. Genov, C. Sun, and X. Zhang, Phys. Rev. Lett., 2008, 85: 123002

[90] D. A. Genov, S. Zhang, and X. Zhang, Nature Physics, 2009, 5: 687

[91] P. Yao, Z. Liang, and X. Jiang, Appl. Phys. Lett., 2008, 92: 031111

[92] H. Chen, Z. Liang, P. Yao, X. Jiang, H. Ma, and C. T. Chan, Phys. Rev. B, 2007, 76: 241104(R)

[93] H. Chen and C. T. Chan, J. Appl. Phys., 2008, 104: 033113

[94] Z. Liang, P. Yao, X. Sun, and X. Jiang, Appl. Phys. Lett., 2008, 92: 131118

[95] B. Zhang, B. I. Wu, and H. Chen, Opt. Express, 2009, 17: 6721

[96] R. Liu, C. Ji, J. J. Mock, J. Y. Chin, T. J. Cui, and D. R. Smith, Science, 2009, 323: 366

[97] J. Valentine, J. Li, T. Zentgraf, G. Bartal, and X. Zhang, Nature Materials, 2009, 8: 568

[98] D. H. Kwon and D. H. Werner, Opt. Express, 2008, 16: 18731

[99] T. Yang, H. Chen, X. Luo, and H. Ma, Opt. Express, 2008, 16: 18545

[100] H. Chen, X. Luo, H. Ma, and C. T. Chan, Opt. Express, 2008, 16: 14603

[101] Y. Lai, H. Chen, Z. Q. Zhang, and C. T. Chan, New J. Phys., 2009, 11: 033010

[102] J. B. Pendry, Phys. Rev. Lett., 2000, 85: 3966

[103] J. B. Pendry and S. A. Ramakrishna, J. Phys.: Condens. Matter, 2002, 14: 8463 8479

[104] J. B. Pendry and S. A. Ramakrishna, J. Phys.: Condens. Matter, 2003, 15: 6345

[105] U. Leonhardt and T. G. Philbin, New J. Phys., 2006, 8: 247

[106] X. Jiang and C. M. Soukoulis, Phys. Rev. Lett., 2000, 85: 70

[107] Y. Lai, J. Ng, H. Chen, D. Z. Han, J. J. Xiao, Z. Q. Zhang, and C. T. Chan, Phys. Rev. Lett., 2009, 102: 253902

Rights & permissions

©Higher Education Press and Springer-Verlag Berlin Heidelberg 2010