Recently, the theory[67, 68] has been developed based on the geometry transformation to realize a cloaking structure (CS), in which objects become invisible from outside. Then a two-dimensional (2D) cylindrical CS[69] and a nonmagnetic optical CS[70, 71] are designed. More surprisingly, the experiment[72] demonstrates that such a 2D 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 invisibility of human dream. However, so far almost all theoretical[67–71, 73–75] 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. [68] and the quantitatively study in our recent work [7], the dispersion is *required* for the cloaking material to avoid the divergent group velocity. For the dispersive CS, new topics, such as the *dynamical process*, can be introduced in. Dynamical study is essential for the cloaking study since without it we can not answer the questions, such as how can the field gets to its stable state?, is there any strong scattering or oscillation in the process?, and how long is the process?, etc. More important, because the real radars generally are pulsive ones, the dynamical process is critical

used to design novel optical/photonic coherence-gain devices.

signal processing.

**dispersive cloaking structures**

Fig. 13. (Color online) (a) The setup of the system. And the distribution of the electric field at different moments during the process. Parameters are chosen as *Ar* = 0.4, *Aθ* = 1.6, *Az* = 0.4 (in the position where *εz* < 1) or *Az* = 1.6 (where *εz* > 1), and *γ* = 0.012*ω*0. (a) *t* = 2.28*T*. (b) *t* = 3.60*T*. (c) *t* = 4.92*T*. (d) *t* = 7.20*T*. (e) *t* = 9.00*T*. (f) Stable state. *T* is the period of the incident EM wave.

for the cloaking effect around the goal frequency. So the dynamical study not only gives us whole physical picture of the cloaking, but also helps us to design more effective cloaks.

In this section, the dynamical process of the electromagnetic (EM) CS is investigated by finite-difference time-domain (FDTD) numerical experiments. In our simulation, the Lorentzian dispersion relations are introduced into the permittivity and the permeability models, then the *real* dynamical process can be simulated[76–78]. Based on numerical simulation, we can follow the details of the dynamical 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 dynamical process and the scattered field, we can find the essential elements which 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. 13(a), similar as the one in Ref. [69]. *R*1 and *R*2 = 2*R*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 *ε*0 = *μ*0 = 1. From the left side, an incident plane wave with working frequency *ω*0 is scattered by the CS, the total field and the scattered field can be recorded inside and outside *B*1 respectively by the numerical technique[79]. So 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 *εz*, *μr*, and *μθ* are needed to be considered (For H-polarized modes, considering the corresponding components *μz*, *εr*, and *εθ*, we can obtain the same numerical results in the dynamical process.). All of them are supposed to have the form 1+ *Fj*(*r*) × *fj*(*ω*), where subscript *j* could be *z*, *r*, and *θ* for *εz*, *μr*, and *μθ*, respectively. The filling factors *Fj*(*r*) are only *r*-dependent, *ωp* is the plasma frequency which set to be a constant *ωp* = 10*ω*0, and *fj*(*ω*) = *ω*2 *p*/(*ω*2 *aj* − *ω*2 − *iωγ*) are the Lorentzian dispersive functions, where

**(d)**

*inct*, (13)

*inct* is the averaged energy flow

**(a) (b)**

**(c)**

The Group Velocity Picture of Metamaterial Systems 79

Fig. 15. (Color online) Direction of Poynting vectors and the intensity front (shown by red dashed curves) at moments during the dynamical process. Parameters are chosen as that of

built up step by step, at last, the field gets to the stable state shown in Fig. 13(f). Because of the dispersion, there is an obvious time delay in the cloaking effect and the strong scattered

We introduce a time-dependent scattering cross-section *σ*(*t*) to quantitatively study the

where *t* = *n* × *T*, *n* = 0, 1, 2, ..., T is the period of the incident wave, ¯*Jscat*(*t*) is the

density of incident field. To observe the dispersive effect on *σ*(*t*) during the dynamical process, at the first step, we keep *Az* and *γ* constant and change *Ar* and *Aθ*, the results are shown in Fig. 14(a). From the *σ* versus *t* curves, we can find the general properties of the dynamical process. First, there is strong scattering in the dynamical process. At the beginning, *σ* increases rapidly when the wave gets 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 dynamical process, called as "relaxation time" generally, can be tuned by the dispersion. From Fig. 14(a), we can see that the main dispersive effect is on the decaying process. From case 1 to case 5, *Ar* and *Aθ* become closer to *one*, so that the dispersion is stronger. We find that the stronger the dispersion, the longer the relaxation time. For comparing with the cloaking cases, we also show the *σ*(*t*) of the naked PEC shell in the case 6. From the definition of *σ*, we know that the area covered by these curve in Fig. 14(a) is proportional to the total scattered energy in the dynamical process. So the CS with *the weaker dispersion will scatter less field* (better cloaking effect) in the dynamical process. But, such

*σ*(*t*) = ¯*Jscat*(*t*)/*S*¯

Fig. 13. (a) *t* = 4.92*T*. (b) *t* = 7.20*T*. (c) *t* = 9.00*T*. (d) Stable state.

one-period-average energy flow of scattered field, and *S*¯

field is observed.

dynamical process, which is defined as

Fig. 14. (Color online) The *σ* vs *t* curves. From (a) to (d), *Az* = 0.4 where *εz* < 1 or *Az* = 1.6 where *εz* > 1. (a) Keep *γ* = 0.012*ω*0 unchanging, choose *Ar* and *Aθ* as case 1: *Ar* = 0.4, *Aθ* = 1.6, case 2: *Ar* = 0.5, *Aθ* = 1.5, case 3: *Ar* = 0.6, *Aθ* = 1.4, case 4: *Ar* = 0.7, *Aθ* = 1.3, case 5: *Ar* = 0.8, *Aθ* = 1.2. Case 6: only PEC shell without CS. (b) Keep *Ar* = 0.4, *Aθ* = 1.6 unchanging, choose *γ*1 = 0.012*ω*0, *γ*2 = 0.024*ω*0, *γ*3 = 0.048*ω*0, *γ*4 = 0.096*ω*0 and *γ*5 = 0.192*ω*0. (c) Keep *Aθ* = 1.6 and *γ* = 0.012*ω*0 unvaried, and change *Ar*. (d) Keep *Ar* = 0.4 and *γ* = 0.012*ω*0 unvaried, and change *Aθ*.

*γ* is the "resonance width" or called as "dissipation factor," *ω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 *εz*, *μr*, and *μθ* always satisfy the geometry transformation of Ref. [69] at *ω*0:

*Re*[*μr*(*r*, *ω*0)] = (*r* − *R*1)/*r*, *Re*[*μθ* (*r*, *ω*0)] = *r*/(*r* − *R*1), and *Re*[*εz*(*r*, *ω*0)] = *R*2 2(*r* − *R*1)/[(*R*2 − *R*1)2*r*]. Then the filling factors *Fj*(*r*) at different *r* can be obtained: *Fr*(*r*) = {*Re*[*μr*(*r*, *ω*0)] − 1}/*Re*[ *fr*(*ω*0)], *Fθ*(*r*) = {*Re*[*μθ*(*r*, *ω*0)] − 1}/*Re*[ *fθ*(*ω*0)], and *Fz*(*r*) = {*Re*[*εz*(*r*, *ω*0)] − 1}/*Re*[ *fz*(*ω*0)].

To investigate the dispersive effect on the dynamical process, we tune the dispersion parameters *ωaj* in our numerical experiments. We use the working frequency *ω*0 as the frequency unit since it is the same for all cases in this section so the ratio *Aj* = *ωaj*/*ω*0 represents *ωaj*. Obviously, for the Lorentzian dispersive relation, the dispersion is stronger when *ω*0 and *ωaj* are closer to each other (the working frequency is near the resonant frequency), or in other words, when *Aj* approaches *one*. Since there are singular values of real part of *ε* and *μ*, in our numerical simulation we have done some approximations,[80] such as we set the maximum and the minimum for *ε* and *μ*. Although such approximations will affect the cloaking effect of stable state,[74] we find that the influence of these approximations on the dynamical process is very small and can be neglected.

First, we show an example of evolving electronic field during the dynamical process in Fig. 13 with concrete parameters of *Ar*, *Aθ*, *Az*, and *γ*. In Fig. 13(a), the plane wave arrives at the left side of the CS, and is ready to enter the CS. From Fig. 13(b)-13(e), the cloaking effect is 22 Will-be-set-by-IN-TECH

**(a)**

**(c)**

**case1 case2 case3 case4 case5 case6**

**Ar =0.4 Ar =0.5 Ar =0.6 Ar =0.7 Ar =0.8** γ **1** γ **2** γ **3** γ **4** γ **5**

**t/T**

Fig. 14. (Color online) The *σ* vs *t* curves. From (a) to (d), *Az* = 0.4 where *εz* < 1 or *Az* = 1.6 where *εz* > 1. (a) Keep *γ* = 0.012*ω*0 unchanging, choose *Ar* and *Aθ* as case 1: *Ar* = 0.4, *Aθ* = 1.6, case 2: *Ar* = 0.5, *Aθ* = 1.5, case 3: *Ar* = 0.6, *Aθ* = 1.4, case 4: *Ar* = 0.7, *Aθ* = 1.3, case 5: *Ar* = 0.8, *Aθ* = 1.2. Case 6: only PEC shell without CS. (b) Keep *Ar* = 0.4, *Aθ* = 1.6

*γ* is the "resonance width" or called as "dissipation factor," *ωaj* are the resonant frequency of

For the study of the dispersive CS, we suppose that the real parts of the *εz*, *μr*, and *μθ* always

*R*1)/[(*R*2 − *R*1)2*r*]. Then the filling factors *Fj*(*r*) at different *r* can be obtained: *Fr*(*r*) = {*Re*[*μr*(*r*, *ω*0)] − 1}/*Re*[ *fr*(*ω*0)], *Fθ*(*r*) = {*Re*[*μθ*(*r*, *ω*0)] − 1}/*Re*[ *fθ*(*ω*0)], and *Fz*(*r*) =

To investigate the dispersive effect on the dynamical process, we tune the dispersion parameters *ωaj* in our numerical experiments. We use the working frequency *ω*0 as the frequency unit since it is the same for all cases in this section so the ratio *Aj* = *ωaj*/*ω*0 represents *ωaj*. Obviously, for the Lorentzian dispersive relation, the dispersion is stronger when *ω*0 and *ωaj* are closer to each other (the working frequency is near the resonant frequency), or in other words, when *Aj* approaches *one*. Since there are singular values of real part of *ε* and *μ*, in our numerical simulation we have done some approximations,[80] such as we set the maximum and the minimum for *ε* and *μ*. Although such approximations will affect the cloaking effect of stable state,[74] we find that the influence of these approximations on

First, we show an example of evolving electronic field during the dynamical process in Fig. 13 with concrete parameters of *Ar*, *Aθ*, *Az*, and *γ*. In Fig. 13(a), the plane wave arrives at the left side of the CS, and is ready to enter the CS. From Fig. 13(b)-13(e), the cloaking effect is

*Re*[*μr*(*r*, *ω*0)] = (*r* − *R*1)/*r*, *Re*[*μθ* (*r*, *ω*0)] = *r*/(*r* − *R*1), and *Re*[*εz*(*r*, *ω*0)] = *R*2

unchanging, choose *γ*1 = 0.012*ω*0, *γ*2 = 0.024*ω*0, *γ*3 = 0.048*ω*0, *γ*4 = 0.096*ω*0 and *γ*5 = 0.192*ω*0. (c) Keep *Aθ* = 1.6 and *γ* = 0.012*ω*0 unvaried, and change *Ar*. (d) Keep

**0 5 10 15 20 25 30**

**A**θ **=1.6 A**θ **=1.5 A**θ **=1.4 A**θ **=1.3 A**θ **=1.2**

**(b)**

**(d)**

2(*r* −

**0 0.5 1.0 1.5 2.0 2.5 3.0 3.5**

**0 0.5 1.0 1.5 2.0 2.5 3.0**

**0 5 10 15 20 25**

*Ar* = 0.4 and *γ* = 0.012*ω*0 unvaried, and change *Aθ*.

satisfy the geometry transformation of Ref. [69] at *ω*0:

the dynamical process is very small and can be neglected.

"atoms" (resonant units) in metamaterials.

{*Re*[*εz*(*r*, *ω*0)] − 1}/*Re*[ *fz*(*ω*0)].

σ**/**λ

Fig. 15. (Color online) Direction of Poynting vectors and the intensity front (shown by red dashed curves) at moments during the dynamical process. Parameters are chosen as that of Fig. 13. (a) *t* = 4.92*T*. (b) *t* = 7.20*T*. (c) *t* = 9.00*T*. (d) Stable state.

built up step by step, at last, the field gets to the stable state shown in Fig. 13(f). Because of the dispersion, there is an obvious time delay in the cloaking effect and the strong scattered field is observed.

We introduce a time-dependent scattering cross-section *σ*(*t*) to quantitatively study the dynamical process, which is defined as

$$
\sigma(t) = \vec{f}\_{\text{scat}}(t) / \vec{S}\_{\text{inct}\prime} \tag{13}
$$

where *t* = *n* × *T*, *n* = 0, 1, 2, ..., T is the period of the incident wave, ¯*Jscat*(*t*) is the one-period-average energy flow of scattered field, and *S*¯ *inct* is the averaged energy flow density of incident field. To observe the dispersive effect on *σ*(*t*) during the dynamical process, at the first step, we keep *Az* and *γ* constant and change *Ar* and *Aθ*, the results are shown in Fig. 14(a). From the *σ* versus *t* curves, we can find the general properties of the dynamical process. First, there is strong scattering in the dynamical process. At the beginning, *σ* increases rapidly when the wave gets 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 dynamical process, called as "relaxation time" generally, can be tuned by the dispersion. From Fig. 14(a), we can see that the main dispersive effect is on the decaying process. From case 1 to case 5, *Ar* and *Aθ* become closer to *one*, so that the dispersion is stronger. We find that the stronger the dispersion, the longer the relaxation time. For comparing with the cloaking cases, we also show the *σ*(*t*) of the naked PEC shell in the case 6. From the definition of *σ*, we know that the area covered by these curve in Fig. 14(a) is proportional to the total scattered energy in the dynamical process. So the CS with *the weaker dispersion will scatter less field* (better cloaking effect) in the dynamical process. But, such

( √ 2*c ε zμr* )/(2 + *ω ε z dε z dω* + *ω μr dμr*

the velocity of light in vacuum.

how to tune *Vgθ* by modifying dispersion parameters.

of the CS, there should be an optimized trade-off.

is performed to validate our analysis.

For the anisotropic cloaking material, the *Vgθ* and *Vgr* can be expressed as: *Vgθ* = [∇**k***ω*(**k**)]*θ* =

The Group Velocity Picture of Metamaterial Systems 81

In order to illustrate our prediction, the *σ*(*t*) under different *Vgθ* and *Vgr* are investigated, respectively. First, we keep the *Vgr* unvaried by holding *Aθ*, *Az* and *γ* constant [keep *dεz*/*dω* and *dμθ*/*dω* unchanged], only modify *Ar* to change the *Vgθ*. The results are shown in Fig.14(c), when *Ar* is closer to *one*, the *Vgθ* becomes smaller (with larger *dμr*/*dω*), the relaxation time is longer and more energy scattered in the dynamical process. So the larger *Vgθ* means the better cloaking effect in the dynamical process. On the other hand, when we keep the *Vgθ* unvaried and change *Vgr* by holding *Ar*, *Az*, and *γ* constant and modifying *Aθ*, the results are shown in Fig. 14(d). We find that the relaxation time is almost unchanged with the change of *Vgr*. Obviously, *Vgθ* is the dominant element in the dynamical process. This conclusion can help us to design a better CS to defend the pulsive radars. In the expression of *Vgθ* , it is also shown

It seems that the larger *Vgθ*, the better cloaking effect in the dynamical process. However, since the *Vg* (and its components) cannot exceed *c* generally, there is a minimum limit for the relaxation time of the cloaking dynamical process. We can estimate it through dividing the mean length of the propagation rays by *Vgθ* . In our model, the mean length is *π*(*R*2 + *R*1)/2, about three wavelengths. So the relaxation time can not be shorter than three periods. Figure 14 shows that our estimation is coincident with our simulation results. Actually, here we are facing a very basic conflict to make a "better" CS, which is more discussed in our other works.[7, 81] The conflict is from the fact that the *pretty strong* dispersion is required to realize a good *stable* cloaking effect at a certain frequency,[7, 68] but at this research we show that the *weaker* dispersion can realize a better cloaking effect in the *dynamical* process. At real design

Based on causality, the limitation of the electromagnetic cloak with dispersive material is investigated in this section The results show that perfect invisibility can not be achieved because of the dilemma that either the group velocity *Vg* 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

Through the ages, people have dreamed to have a magic cloak whose owner can not be seen by others. For this fantastic dream, plenty of work has been done by scientists all over the world. For example, the researchers diminished the scattering or the reflection from objects by absorbing screens[82] and small, non-absorbing, compound ellipsoids[83]. More recently, based on the coordinate transformation, J. B. Pendry, et al theoretically proposed a general recipe for designing an electromagnetic cloak to hide an object from the electromagnetic(EM) wave[68]. 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[67, 75]. Soon, Steven A. Cummer, et al simulated numerically(COMSOL) the cylindrical version of this cloak structure using ideal and nonideal (but physically realizable) electromagnetic

**6. Limitation of the electromagnetic cloak with dispersive material**

2*c ε zμθ* )/(2 + *ω ε z dε z dω* + *ω μθ dμθ*

*dω* ), where *c* is

*dω* ) and *Vgr* = [∇**k***ω*(**k**)]*r* = ( √

a general conclusion is still not enough for us to get a clear physical picture to understand the cloaking dynamical process.

Next, we check whether the absorption of the CS is important in the process. The absorption is determined by the imaginary part of *ε* and *μ*. To study this effect, we hold *Ar*, *Aθ*, and *Az* constant but modify the dispassion factor *γ*. We modify the filling factors *Fj* simultaneously, so that the real parts of *ε* and *μ* are kept unchanged at *ω*0. In such way, we can keep the dispersion strength almost unchanged, but with the imaginary parts of *ε* and *μ* changed. Results in Fig. 14(b) show that the stronger absorption only leads 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 *σ*(*t*) of stable state considerably.

To obtain deeper insight of the dynamical process, we need to study the dynamical process more carefully. From Figs. 13(b)-13(e), we can see that the "field intensity" (shown by different color in the figures) propagates slower inside the CS than that in the outside vacuum. And when the inside field intensity "catches up" the outside one [in Fig. 13(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 dynamical process. From Figs. 15(a)-15(d), we show the direction of Poynting vectors in moments of Figs. 13(c)-13(f), respectively. In the Fig. 15, we see that there is the "intensity front"(shown by red dashed curve) which separates 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). But at the left-side region which is swept by the intensity front, the Poynting vectors are very regular and nearly along the "cloaking rays" which was predicted at the coordinate transformation[68]. Since 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 is related the original cloaking recipe,[68] which makes the cloaking material is almost impedance matched layer by layer. This also explains why there is no oscillation in the cloaking dynamical process generally. This picture also can interpret the strong scattered field in the dynamical process, since these "irregular rays" at the right-side region of the intensity front must be scattered strongly. Further, we can use this picture to analysis the dynamical process of other incident waves, such as the Gaussian beams, which are composed of different plane-wave components.

With these 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 *Vg*, which is controlled by the material dispersion. So the intensity front, which determines the dynamical process, should move in *Vg*. Thus, we can explain the results in Fig. 14, since our modification of the dispersive parameters can cause the *Vg* changed. But, because the cloaking material is the strong anisotropic material, the *Vg* at different directions could be very different. Can we predict more precisely which component dominates the relaxation time? The answer is "yes." In Fig. 15(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 *Vgθ*, not the component along the *r* direction *Vgr*, that dominates the relaxation time and the total scattered energy in the dynamical process.

24 Will-be-set-by-IN-TECH

a general conclusion is still not enough for us to get a clear physical picture to understand the

Next, we check whether the absorption of the CS is important in the process. The absorption is determined by the imaginary part of *ε* and *μ*. To study this effect, we hold *Ar*, *Aθ*, and *Az* constant but modify the dispassion factor *γ*. We modify the filling factors *Fj* simultaneously, so that the real parts of *ε* and *μ* are kept unchanged at *ω*0. In such way, we can keep the dispersion strength almost unchanged, but with the imaginary parts of *ε* and *μ* changed. Results in Fig. 14(b) show that the stronger absorption only leads to larger stable value of *σ*, leaving the relaxation time nearly unchanged. Thus, we can exclude the absorption from the relevant

To obtain deeper insight of the dynamical process, we need to study the dynamical process more carefully. From Figs. 13(b)-13(e), we can see that the "field intensity" (shown by different color in the figures) propagates slower inside the CS than that in the outside vacuum. And when the inside field intensity "catches up" the outside one [in Fig. 13(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 dynamical process. From Figs. 15(a)-15(d), we show the direction of Poynting vectors in moments of Figs. 13(c)-13(f), respectively. In the Fig. 15, we see that there is the "intensity front"(shown by red dashed curve) which separates 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). But at the left-side region which is swept by the intensity front, the Poynting vectors are very regular and nearly along the "cloaking rays" which was predicted at the coordinate transformation[68]. Since 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 is related the original cloaking recipe,[68] which makes the cloaking material is almost impedance matched layer by layer. This also explains why there is no oscillation in the cloaking dynamical process generally. This picture also can interpret the strong scattered field in the dynamical process, since these "irregular rays" at the right-side region of the intensity front must be scattered strongly. Further, we can use this picture to analysis the dynamical process of other incident waves, such as the Gaussian beams, which are composed of different

With these 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 *Vg*, which is controlled by the material dispersion. So the intensity front, which determines the dynamical process, should move in *Vg*. Thus, we can explain the results in Fig. 14, since our modification of the dispersive parameters can cause the *Vg* changed. But, because the cloaking material is the strong anisotropic material, the *Vg* at different directions could be very different. Can we predict more precisely which component dominates the relaxation time? The answer is "yes." In Fig. 15(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 *Vgθ*, not the component along the *r* direction *Vgr*, that dominates the relaxation time and the total scattered energy in the dynamical process.

parameter list, since it only influences the *σ*(*t*) of stable state considerably.

cloaking dynamical process.

plane-wave components.

For the anisotropic cloaking material, the *Vgθ* and *Vgr* can be expressed as: *Vgθ* = [∇**k***ω*(**k**)]*θ* = ( √ 2*c ε zμr* )/(2 + *ω ε z dε z dω* + *ω μr dμr dω* ) and *Vgr* = [∇**k***ω*(**k**)]*r* = ( √ 2*c ε zμθ* )/(2 + *ω ε z dε z dω* + *ω μθ dμθ dω* ), where *c* is the velocity of light in vacuum.

In order to illustrate our prediction, the *σ*(*t*) under different *Vgθ* and *Vgr* are investigated, respectively. First, we keep the *Vgr* unvaried by holding *Aθ*, *Az* and *γ* constant [keep *dεz*/*dω* and *dμθ*/*dω* unchanged], only modify *Ar* to change the *Vgθ*. The results are shown in Fig.14(c), when *Ar* is closer to *one*, the *Vgθ* becomes smaller (with larger *dμr*/*dω*), the relaxation time is longer and more energy scattered in the dynamical process. So the larger *Vgθ* means the better cloaking effect in the dynamical process. On the other hand, when we keep the *Vgθ* unvaried and change *Vgr* by holding *Ar*, *Az*, and *γ* constant and modifying *Aθ*, the results are shown in Fig. 14(d). We find that the relaxation time is almost unchanged with the change of *Vgr*. Obviously, *Vgθ* is the dominant element in the dynamical process. This conclusion can help us to design a better CS to defend the pulsive radars. In the expression of *Vgθ* , it is also shown how to tune *Vgθ* by modifying dispersion parameters.

It seems that the larger *Vgθ*, the better cloaking effect in the dynamical process. However, since the *Vg* (and its components) cannot exceed *c* generally, there is a minimum limit for the relaxation time of the cloaking dynamical process. We can estimate it through dividing the mean length of the propagation rays by *Vgθ* . In our model, the mean length is *π*(*R*2 + *R*1)/2, about three wavelengths. So the relaxation time can not be shorter than three periods. Figure 14 shows that our estimation is coincident with our simulation results. Actually, here we are facing a very basic conflict to make a "better" CS, which is more discussed in our other works.[7, 81] The conflict is from the fact that the *pretty strong* dispersion is required to realize a good *stable* cloaking effect at a certain frequency,[7, 68] but at this research we show that the *weaker* dispersion can realize a better cloaking effect in the *dynamical* process. At 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 can not be achieved because of the dilemma that either the group velocity *Vg* 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.

#### **6. Limitation of the electromagnetic cloak with dispersive material**

Through the ages, people have dreamed to have a magic cloak whose owner can not be seen by others. For this fantastic dream, plenty of work has been done by scientists all over the world. For example, the researchers diminished the scattering or the reflection from objects by absorbing screens[82] and small, non-absorbing, compound ellipsoids[83]. More recently, based on the coordinate transformation, J. B. Pendry, et al theoretically proposed a general recipe for designing an electromagnetic cloak to hide an object from the electromagnetic(EM) wave[68]. 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[67, 75]. Soon, Steven A. Cummer, et al simulated numerically(COMSOL) the cylindrical version of this cloak structure using ideal and nonideal (but physically realizable) electromagnetic

Fig. 17. The relation between *R*1/*r* and *α* when *Vg* = *c* for different *ω*

the dispersion relation of the anisotropic material[87] as: *k*2

(*cosα*)2 *n*2 *r*�

*dω* and *mt*� = *nt*� + *ω*

We will discuss Eq.15 in two cases. The first case is with the finite *dnt*�

tangential component of *Vg*) will diverge at both peaks around *α* = 0.

+ (*sinα*)2 *n*2 *t*�

*dnt*� *dω* . If the transformation has the following characteristic: *f*(*r* = 0) = *R*1, *f*(*r* = *R*2) = *R*2, then when *r*� → *R*1 (or *r* → 0), *nt*� will tend to zero, and the group velocity is approximated as:

*Vg* ≈ *c*

example, in which the transformation is *r*� = *f*(*r*)=(*R*2 − *R*1)*r*/*R*2 + *R*1 as Ref. [68], *R*2 = 2*R*1, thus *nr*� = 2 and *nt*� = 2 − 4/(*r*/*R*1 + 2). The dispersion parameters are set as *mr*� = 2.5,

*dω* = 4 at working frequency. In figure 16, the curves of *Vg* vs *α* are plotted for different *R*1/*r* values. We can see that, for large *R*1/*r*(*r* → 0), the group velocity (more precisely, the

*ϕ*� , *nr*� = √*εϕ*�*μθ*� = *ni*/ *d f*(*r*)

*k*2 *t*� = *k*2

*ω dnt*� *θ*� + *k*2

Where *mr*� = *nr*� + *ω*

obtained as:

*ni* = √*εiμi*. Then we can define *kr*� = *ω*

*Vg* = *c*

*dnr*�

diverge when *sinα* → 0 for any finite *dnt*�

is very important for us to have consistent physical pictures in both spaces. At working frequency *ω*0, for a propagating mode with k-vector as {*kr*� ,*kθ*� ,*kϕ*�} inside the cloak, we have

The Group Velocity Picture of Metamaterial Systems 83

*c nr*� *cosα* and *kt*� = *ω*

/((*cosα*)2 *mr*�

*dnt*� *dω*

*nr*�

*dnt*� *dω* .

*c nt*�*sinα*, the group velocity can be

*t*� = *ω*2/*c*2, where

) (14)

*dω* . Obviously, *Vg* will

(15)

*r*�/*n*2 *r*� + *k*2 *t*�/*n*2

*dr* , *nt*� = *nθ*� = *nϕ*� = √*εr*�*μθ*� = *nir*/ *f*(*r*), and

+ (*sinα*)2 *mt*�

*dω* . Such divergence is shown in figure 16 for a concrete

*nt*�

Fig. 16. The group velocity *Vg* versus *α* for different *R*1/*r* values

parameters[69]. Especially, Schurig, et al experimentally demonstrated such a cloak by split-ring resonators[72]. In addition, Wenshan 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[70, 71]. According to the general recipe, the electromagnetic cloak is supposed to be perfect or "fully functioned" at 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 [84, 85]. 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 very basic physical laws, are not well studied. If the dispersion is introduced into the study, can we have a deeper insight into the cloaking physics?

In this section we will show the ideal cloaking can not be achieved because of another more basic physical limitation—the causality limitation (based on the same limitation, Chen et al obtained a constraint of the band width that limit the design of an invisibility cloak[86]). Starting from dispersion relation and combining with the demand of causality, we will demonstrate that the ideal cloaking will lead to the dilemma that either the group velocity *Vg* diverges or a strong absorption is imposed on the cloaking material. Our derivation and numerical experiments based on the finite-difference-time-domain(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's consider a more general coordinate transformation on an initial homogeneous medium with *εi* = *μi* in r space: *r* = *f*(*r*), *θ* = *θ*, *ϕ* = *ϕ*, following the approach in Ref. [87] and [89], we get the following radius-dependent, anisotropic relative permittivity and permeability: *εr* = *μr* = *εi*( *r f*(*r*))2 *d f*(*r*) *dr* , *εθ* = *μθ* = *εi*/ *d f*(*r*) *dr* and *εϕ* = *μϕ* = *εi*/ *d f*(*r*) *dr* . 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 physical properties of wave propagation in *r* space should be inherited in *r* space, such as the absorption. This 26 Will-be-set-by-IN-TECH

R1 /r=0.5

R1 /r=1.5

R1 /r=3.5

R1 /r=7.5

R1 /r=15.5

R1 /r=18.5

R1 /r=31.5

R1 /r=63.5

−1.5 −1 −0.5 0 0.5 1 1.5

parameters[69]. Especially, Schurig, et al experimentally demonstrated such a cloak by split-ring resonators[72]. In addition, Wenshan 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[70, 71]. According to the general recipe, the electromagnetic cloak is supposed to be perfect or "fully functioned" at 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 [84, 85]. 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 very basic physical laws, are not well studied. If the dispersion is introduced into the study, can we have a deeper

In this section we will show the ideal cloaking can not be achieved because of another more basic physical limitation—the causality limitation (based on the same limitation, Chen et al obtained a constraint of the band width that limit the design of an invisibility cloak[86]). Starting from dispersion relation and combining with the demand of causality, we will demonstrate that the ideal cloaking will lead to the dilemma that either the group velocity *Vg* diverges or a strong absorption is imposed on the cloaking material. Our derivation and numerical experiments based on the finite-difference-time-domain(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's consider a more general coordinate transformation on an initial homogeneous medium with *εi* = *μi* in r space: *r* = *f*(*r*), *θ* = *θ*, *ϕ* = *ϕ*, following the approach in Ref. [87] and [89], we get the following radius-dependent, anisotropic relative permittivity and permeability:

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 physical properties of wave propagation in *r* space should be inherited in *r* space, such as the absorption. This

*dr* and *εϕ* = *μϕ* = *εi*/ *d f*(*r*)

*dr* . We emphasize that

α/rad

0.2

insight into the cloaking physics?

*f*(*r*))2 *d f*(*r*)

*dr* , *εθ* = *μθ* = *εi*/ *d f*(*r*)

*εr* = *μr* = *εi*( *r*

Fig. 16. The group velocity *Vg* versus *α* for different *R*1/*r* values

0.6

1

Vg/c

1.4

Fig. 17. The relation between *R*1/*r* and *α* when *Vg* = *c* for different *ω dnt*� *dω* .

is very important for us to have consistent physical pictures in both spaces. At working frequency *ω*0, for a propagating mode with k-vector as {*kr*� ,*kθ*� ,*kϕ*�} inside the cloak, we have the dispersion relation of the anisotropic material[87] as: *k*2 *r*�/*n*2 *r*� + *k*2 *t*�/*n*2 *t*� = *ω*2/*c*2, where *k*2 *t*� = *k*2 *θ*� + *k*2 *ϕ*� , *nr*� = √*εϕ*�*μθ*� = *ni*/ *d f*(*r*) *dr* , *nt*� = *nθ*� = *nϕ*� = √*εr*�*μθ*� = *nir*/ *f*(*r*), and *ni* = √*εiμi*. Then we can define *kr*� = *ω c nr*� *cosα* and *kt*� = *ω c nt*�*sinα*, the group velocity can be obtained as:

$$V\_{\mathcal{S}} = c\sqrt{\frac{(\cos a)^2}{n\_{r'}^2} + \frac{(\sin a)^2}{n\_{l'}^2}}/((\cos a)^2 \frac{m\_{r'}}{n\_{r'}} + (\sin a)^2 \frac{m\_{l'}}{n\_{l'}})\tag{14}$$

Where *mr*� = *nr*� + *ω dnr*� *dω* and *mt*� = *nt*� + *ω dnt*� *dω* .

If the transformation has the following characteristic: *f*(*r* = 0) = *R*1, *f*(*r* = *R*2) = *R*2, then when *r*� → *R*1 (or *r* → 0), *nt*� will tend to zero, and the group velocity is approximated as:

$$V\_{\mathcal{S}} \approx \frac{c}{|\sin \mathfrak{a}| \omega \frac{d n\_{\ell'}}{d \omega}} \tag{15}$$

We will discuss Eq.15 in two cases. The first case is with the finite *dnt*� *dω* . Obviously, *Vg* will diverge when *sinα* → 0 for any finite *dnt*� *dω* . Such divergence is shown in figure 16 for a concrete example, in which the transformation is *r*� = *f*(*r*)=(*R*2 − *R*1)*r*/*R*2 + *R*1 as Ref. [68], *R*2 = 2*R*1, thus *nr*� = 2 and *nt*� = 2 − 4/(*r*/*R*1 + 2). The dispersion parameters are set as *mr*� = 2.5, *ω dnt*� *dω* = 4 at working frequency. In figure 16, the curves of *Vg* vs *α* are plotted for different *R*1/*r* values. We can see that, for large *R*1/*r*(*r* → 0), the group velocity (more precisely, the tangential component of *Vg*) will diverge at both peaks around *α* = 0.

(a) (b)

structure with a PEC at radius *R*1, (b) the naked PEC with radius *R*1.

needed to avoid *Vg* > *c*.

*�*˜*z*�(*r*�

*μ*˜*r*�(*r*�

*μ*˜ *θ*�(*r*�

, *ω*)| = 1 + *Fz*�(*r*�

, *ω*) = 1 + *Fr*�(*r*�

, *ω*) = 1 + *Fθ*�(*r*�

Fig. 18. The snapshots of the electric-field distribution in the vicinity of PEC. (a) the cloaking

The Group Velocity Picture of Metamaterial Systems 85

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. [68], the cloaking material (in *r*� space) must be dispersive to avoid the group velocity over c. Such Lorentz co-variant violation is generally true for ˛aˇrtransformation optics ˛a´s since material parameters are non-relativistic, so the causality limitation should be checked widely. Third, from Eq.14, we can find that not only the inner layers of the cloak(*r*� → *R*1) but also the other layers(*r*� > *R*1) must be dispersive. For every layer, a certain dispersive strength is

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 two-dimensional cloak[69]. Without lost of generality, only TE modes are investigated in this study (TE modes have the electric field perpendicular to the two-dimensional plane of our model). Thus the constitutive parameters involved here are *�z*� , *μr*� , *μθ*� . The dispersion is introduced into our FDTD by standard Lorentz model:

)*ω*2

)*ω*2

)*ω*2

Where *ωpz*� , *ωpr*� , *ωpθ*� are plasma frequencies, *ωaz*� , *ωar*�, *ωaθ*� are atom resonated frequencies, *γz*� , *γr*� , *γθ*� are damping factors and *Fpz*� , *Fpr*� , *Fpθ*� are filling factors. In our simulation, an E-polarized time-harmonic uniform plane wave whose wavelength *λ*0 in vacuum is 3.75*cm* is incident from left to right. The real parts of the constitutive parameters at *ω*0 = 2*πc*/*λ*0

*r*� = *R*2(*θ*) − *R*1(*θ*)

*pz*�/(*ωaz*�(*r*�

*pr*�/(*ωar*�(*r*�

*pθ*�/(*ωaθ*�(*r*�

)

)

)

2 − *ω*2 − *iωγz*�)

2 − *ω*2 − *iωγθ*�)

*θ*� = *θ* (18)

*R*2(*θ*) *r* + *R*1(*θ*) (17)

2 − *ω*2 − *iωγr*�) (16)

Because of the causality limitation, it is well-known that the group velocity can not exceed *c* except in the "strong dispersion" frequency range (or called "resonant range"). But, 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. So perfect invisibility can not be achieved for the finite *dnt*� *dω* because it will lead to superluminal velocity or strong absorption.

In addition, the curves with the criterion condition *Vg* = *c* on the plane [*R*1/*r*, *α*] are plotted for different *ω dnt*� *dω* in figure 17. The region to the left of curves is corresponding to *Vg* < *c* and the region to the right is corresponding to *Vg* > *c*. There exists a maximum *max*{*R*1/*r*} for each curve in order that *Vg* ≤ *c* can be hold for all the *α*. Especially, for the no-dispersion case *ω dnt*� *dω* = 0, we can see that *Vg* > *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.14, and it is consistent with the analysis in Ref. [68]. From figure 17, we know that the larger *ω dnt*� *dω* , the larger *max*{*R*1/*r*}. But anyway, for arbitrary finite *ω dnt*� *dω* , *max*{*R*1/*r*} can not be infinite, so that the superluminal range always exists.

The second case of Eq.15 is with divergent *dnt*� *dω* . From the previous discussion, we know that if the ideal cloak exists, the cloak must be dispersive and *dnt*� *dω* must be divergent when *r* → 0. Actually, when *r* → 0, since √*εr*� ∝ *r*, *dnt*� *dω* ∝ *dεr*� *dω* / √*εr*� is really divergent for non-zero *dεr*� *dω* . From Eq.15, we can see that now the *Vg* → 0 for a finite *dεr*� *dω* (generally true) at all *α* values except *α* = 0(or *π*), so that the group velocity difficulty seems to be overcome. But, since the causality limitation, the non-zero *dεr*� *dω* 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*1 is almost totally absorbed by the cloaking material. As we pointed out at the beginning that the absorption in *r*� space should also appear in r space, because of the consistence 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* → 0, the finite imaginary part in *r*� space corresponds to the infinite imaginary part in *r* space, which also means very strong absorption in the initial homogeneous medium. So the perfect cloaking is still impossible because of the strong absorption which is enforced by the causality limitation.

For a two-dimensional coordinate transformation: *r*� = *f*(*r*), *θ*� = *θ*, *z*� = *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 *Vg* diverges or a strong absorption is imposed on the cloaking material. First, it is an intrinsic conflict which can not be solved by the methods, for example, "the system is imbedded in a medium"[68]. We believe that the ideal cloaking is impossible because of the causality limitation and this conclusion is consistent with the statement of previous studies[89] that the perfect invisibility is unachievable because of the wave nature of light. Second, we have to face the question: " Why the causality is violated for ideal cloaking which is based on the simple coordinate transformation?". Our answer is that the causality is only guaranteed by 28 Will-be-set-by-IN-TECH

Because of the causality limitation, it is well-known that the group velocity can not exceed *c* except in the "strong dispersion" frequency range (or called "resonant range"). But, 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. So perfect invisibility

In addition, the curves with the criterion condition *Vg* = *c* on the plane [*R*1/*r*, *α*] are plotted

and the region to the right is corresponding to *Vg* > *c*. There exists a maximum *max*{*R*1/*r*} for each curve in order that *Vg* ≤ *c* can be hold for all the *α*. Especially, for the no-dispersion

*dω* . From Eq.15, we can see that now the *Vg* → 0 for a finite *dεr*�

true) at all *α* values except *α* = 0(or *π*), so that the group velocity difficulty seems to be

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*1 is almost totally absorbed by the cloaking material. As we pointed out at the beginning that the absorption in *r*� space should also appear in r space, because of the consistence 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* → 0, the finite imaginary part in *r*� space corresponds to the infinite imaginary part in *r* space, which also means very strong absorption in the initial homogeneous medium. So the perfect cloaking is still impossible because of the strong absorption which is enforced

For a two-dimensional coordinate transformation: *r*� = *f*(*r*), *θ*� = *θ*, *z*� = *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

Next, we will discuss the physical meaning of the dilemma that either the group velocity *Vg* diverges or a strong absorption is imposed on the cloaking material. First, it is an intrinsic conflict which can not be solved by the methods, for example, "the system is imbedded in a medium"[68]. We believe that the ideal cloaking is impossible because of the causality limitation and this conclusion is consistent with the statement of previous studies[89] that the perfect invisibility is unachievable because of the wave nature of light. Second, we have to face the question: " Why the causality is violated for ideal cloaking which is based on the simple coordinate transformation?". Our answer is that the causality is only guaranteed by

*dω* = 0, we can see that *Vg* > *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.14, and it is consistent with the analysis in Ref. [68]. From figure 17,

*dω* in figure 17. The region to the left of curves is corresponding to *Vg* < *c*

*dω* , the larger *max*{*R*1/*r*}. But anyway, for arbitrary finite *ω*

*dω* ∝ *dεr*�

*dω* /

*dω* . From the previous discussion, we know

*dω* must be divergent when

√*εr*� is really divergent for

*dω* means non-zero imaginary part

*dω* (generally

*dω* because it will lead to superluminal velocity or strong

*dnt*� *dω* ,

can not be achieved for the finite *dnt*�

*dnt*�

absorption.

case *ω*

non-zero

case.

for different *ω*

*dnt*�

we know that the larger *ω*

*dεr*�

by the causality limitation.

*dnt*�

*r* → 0. Actually, when *r* → 0, since √*εr*� ∝ *r*, *dnt*�

overcome. But, since the causality limitation, the non-zero *dεr*�

The second case of Eq.15 is with divergent *dnt*�

*max*{*R*1/*r*} can not be infinite, so that the superluminal range always exists.

that if the ideal cloak exists, the cloak must be dispersive and *dnt*�

Fig. 18. The snapshots of the electric-field distribution in the vicinity of PEC. (a) the cloaking structure with a PEC at radius *R*1, (b) the naked PEC with radius *R*1.

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. [68], the cloaking material (in *r*� space) must be dispersive to avoid the group velocity over c. Such Lorentz co-variant violation is generally true for ˛aˇrtransformation optics ˛a´s since material parameters are non-relativistic, so the causality limitation should be checked widely. Third, from Eq.14, we can find that not only the inner layers of the cloak(*r*� → *R*1) but also the other layers(*r*� > *R*1) must be dispersive. For every layer, a certain dispersive strength is needed to avoid *Vg* > *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 two-dimensional cloak[69]. Without lost of generality, only TE modes are investigated in this study (TE modes have the electric field perpendicular to the two-dimensional plane of our model). Thus the constitutive parameters involved here are *�z*� , *μr*� , *μθ*� . The dispersion is introduced into our FDTD by standard Lorentz model:

$$\begin{aligned} \left| \mathfrak{E}\_{z'} (r', \omega) \right| &= 1 + F\_{z'} (r') \omega\_{p z'}^2 / (\omega\_{a z'} (r')^2 - \omega^2 - \mathrm{i}\omega \gamma\_{z'})\\ \left| \mathfrak{H}\_{r'} (r', \omega) = 1 + F\_{r'} (r') \omega\_{p r'}^2 / (\omega\_{a r'} (r')^2 - \omega^2 - \mathrm{i}\omega \gamma\_{r'})\\ \left| \mathfrak{H}\_{\theta'} (r', \omega) = 1 + F\_{\theta'} (r') \omega\_{p \theta'}^2 / (\omega\_{a \theta'} (r')^2 - \omega^2 - \mathrm{i}\omega \gamma\_{\theta'}) \end{aligned} \tag{16}$$

$$r' = \frac{R\_2(\theta) - R\_1(\theta)}{R\_2(\theta)}r + R\_1(\theta) \tag{17}$$

$$
\theta' = \theta \tag{18}
$$

Where *ωpz*� , *ωpr*� , *ωpθ*� are plasma frequencies, *ωaz*� , *ωar*�, *ωaθ*� are atom resonated frequencies, *γz*� , *γr*� , *γθ*� are damping factors and *Fpz*� , *Fpr*� , *Fpθ*� are filling factors. In our simulation, an E-polarized time-harmonic uniform plane wave whose wavelength *λ*0 in vacuum is 3.75*cm* is incident from left to right. The real parts of the constitutive parameters at *ω*0 = 2*πc*/*λ*0

metamaterial studies. From group velocity, we can find the physical origin of abnormal optical phenomena of metamaterials, such as the "no-transmission no-reflection" on the hyper-medium surface which is from a zero-group-velocity reflecting mode, and the coherence gain of superlens image which is from the different group delay on different paths. From group velocity, we can avoid some traps of violating basic physical limitation, such as the violation of causality limitation in cloaking study. These traps are very serious since the metamaterial from our imagination could exist in this world if violating basic limitations. From group velocity, we can find the key parameter of cloaking dynamical process and help us to optimizing the design of cloak design. From group velocity of evanescent wave, new detecting methods could be found, for example the detecting of QED vacuum polarization by phase change or delay time of evanescent wave. We believe that only with the well-constructed group velocity picture, the deeper understanding of the abnormal optical/photonic properties of metamaterials is possible. All these research works also show that the group velocity study of metamaterials can lead us to many new interesting topics,

The Group Velocity Picture of Metamaterial Systems 87

This work is supported by NKBRPC (Grant No. 2012CB927401), NSFC (Grant Nos. 11004212, 11174309, and 60938004), and the STCSM (Grant Nos. 11ZR1443800 and 11JC1414500).

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

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

[5] J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, IEEE Trans. Microwave

[6] P. Markos and C. M. Soukoulis, Phys. Rev. B 65, 033401 (2001); Phys. Rev. E 65, 036622

[8] Zixian Liang, Peijun Yao, Xiaowei Sun and Xunya Jiang, Appl. Phys. Lett. 92, 131118

[9] Zheng Liu, Zixian Liang, Xunya Jiang, Xinhua Hu, Xin Li, and Jian Zi, Appl.Phys. Lett.

[10] Brillouin, *Wave Propagation and Group Velocity*, New York:Academic Press Inc., (1960) [11] J. H. Page, Ping Sheng, H. P. Schriemer, I.Jones, Xiaodun Jing and D. A. Weitzt,Science,

[13] R.A. Shelby, D.R. Smith, and S. Schultz, Science 92, 792(2001); D.R. Smith and D. Schurig, Phys. Rev. Lett. 90,077405 (2003);C.G. Parazzoli, R.B. Greegor, K. Li, B.E.C.

[14] J. B. Pendry, Phys. Rev. Lett. 85, 3966(2000); N. Fang, H. Lee, C. Sun, X. Zhang, Science 308, 534 (2005); T. Taubner, D. Korobkin, Y. Urzhumov, G. Shvets, R.Hillenbrand, Science

Theory Tech. 47, 2075 (1999); J. B. Pendry, Phys. Rev. Lett. 85, 3966 (2000).

which are still waiting for further research.

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

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

[7] P. Yao, Z. Liang, and X. Jiang, Appl. Phys. Lett. 92, 031111 (2008).

[12] Hailu Luo, Wei Hu and Zhongzhou Ren, Europhys. Lett. 74 1081(2006).

Koltenbah, and M. Tanielian, Phys. Rev. Lett. 90, 107401 (2003).

**8. Acknowledgement**

84, 4184 (2000).

**9. References**

(2001).

(2002).

(2008)

96, 113507 (2010)

271 634 1996

313, 1595 (2006)

satisfy the cloaking coordinate transformation[69, 91], and they are *μr*� = *r*� −*R*1 *r*� , *μθ*� = 1 *μr*� , *εz*� = ( *R*2 *R*2−*R*1 )2 *r*� −*R*1 *r*� , where *R*1 is 0.665*λ*0, *R*2 is 1.33*λ*0. And the dispersive parameters are set as follows: if *εz*� > 1 then *ωaz*� = 1.4*ω*0, else *ωaz*� = 0.6*ω*0, *ωar*� = 0.6*ω*0, *ωaθ*� = 1.4*ω*0, *γz*� = *γr*� = *γθ*� = *ω*0/100. *ωpz*� = *ωpr*� = *ωpθ*� = 4*ω*0, *Fz*�(*r*)=(*εz*� − 1) (*ω*2 *az*�−*ω*2 0 )2+*ω*2 0*γ*2 *z*� (*ω*2 *az*�−*ω*2 0 )*ω*2 *pz*� , *Fr*�(*r*) = (*μr*� − 1) (*ω*2 *ar*�−*ω*2 0 )2+*ω*2 0*γ*2 *r*� (*ω*2 *ar*�−*ω*2 0 )*ω*2 *pr*� , *Fθ*�(*r*)=(*μθ*� − 1) (*ω*2 *aθ*�−*ω*2 0 )2+*ω*2 0*γ*2 *θ*� (*ω*2 *aθ*�−*ω*2 0 )*ω*2 *pθ*� . 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 [92].

Figure 18 shows the snapshots of the electric-field distribution in two cases: the cloak with the perfect electric conductor (PEC) at radius *R*1(left), and the naked PEC with radius *R*1(right). 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*λ*0 and 0.24*λ*0 respectively, while the scattering cross section of the naked PEC is 3.14*λ*0. So, 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 *ε* = *μ* = 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*λ*0 which is about one order smaller. The reason of strong absorption has been discussed before.

Now we can have a full view of cloaking recipe based on the coordinate transformation. First, the cloaking material must be dispersive, and the strong absorption can not 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 that the scattered field is weak. Although the ideal invisibility is impossible, the cloaking recipe still has a main advantage. The "strong-absorption and weak-scattering" property means that the cloak almost can not be observed except from the forward direction. so that such cloak can well defend the aˇ˛rpassive radars ˛a´s 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 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 a strong absorption can not be overcome since 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.

#### **7. Summary**

In summary, we have investigate the metamaterial systems from group velocity(energy velocity) picture. From these topics, we demonstrate the importance of group velocity in metamaterial studies. From group velocity, we can find the physical origin of abnormal optical phenomena of metamaterials, such as the "no-transmission no-reflection" on the hyper-medium surface which is from a zero-group-velocity reflecting mode, and the coherence gain of superlens image which is from the different group delay on different paths. From group velocity, we can avoid some traps of violating basic physical limitation, such as the violation of causality limitation in cloaking study. These traps are very serious since the metamaterial from our imagination could exist in this world if violating basic limitations. From group velocity, we can find the key parameter of cloaking dynamical process and help us to optimizing the design of cloak design. From group velocity of evanescent wave, new detecting methods could be found, for example the detecting of QED vacuum polarization by phase change or delay time of evanescent wave. We believe that only with the well-constructed group velocity picture, the deeper understanding of the abnormal optical/photonic properties of metamaterials is possible. All these research works also show that the group velocity study of metamaterials can lead us to many new interesting topics, which are still waiting for further research.

#### **8. Acknowledgement**

This work is supported by NKBRPC (Grant No. 2012CB927401), NSFC (Grant Nos. 11004212, 11174309, and 60938004), and the STCSM (Grant Nos. 11ZR1443800 and 11JC1414500).

#### **9. References**

30 Will-be-set-by-IN-TECH

as follows: if *εz*� > 1 then *ωaz*� = 1.4*ω*0, else *ωaz*� = 0.6*ω*0, *ωar*� = 0.6*ω*0, *ωaθ*� = 1.4*ω*0, *γz*� =

(*ω*2 *aθ*�−*ω*2

many possible choices. The different groups of parameters correspond to different dynamic

Figure 18 shows the snapshots of the electric-field distribution in two cases: the cloak with the perfect electric conductor (PEC) at radius *R*1(left), and the naked PEC with radius *R*1(right). 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*λ*0 and 0.24*λ*0 respectively, while the scattering cross section of the naked PEC is 3.14*λ*0. So, 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 *ε* = *μ* = 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*λ*0 which is about one order smaller. The reason of strong absorption has been discussed

Now we can have a full view of cloaking recipe based on the coordinate transformation. First, the cloaking material must be dispersive, and the strong absorption can not 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 that the scattered field is weak. Although the ideal invisibility is impossible, the cloaking recipe still has a main advantage. The "strong-absorption and weak-scattering" property means that the cloak almost can not be observed except from the forward direction. so that such cloak can well defend the aˇ˛rpassive radars ˛a´s 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 giant Rayleigh scattering case,

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 a strong absorption can not be overcome since 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

In summary, we have investigate the metamaterial systems from group velocity(energy velocity) picture. From these topics, we demonstrate the importance of group velocity in

(*ω*2 *aθ*�−*ω*2 0 )*ω*2 *pθ*�

*r*� , where *R*1 is 0.665*λ*0, *R*2 is 1.33*λ*0. And the dispersive parameters are set

0 )2+*ω*2 0*γ*2 *θ*� −*R*1

0 )2+*ω*2 0*γ*2 *z*�

. In fact, these parameters have

(*ω*2 *az*�−*ω*2

(*ω*2 *az*�−*ω*2 0 )*ω*2 *pz*�

*r*� , *μθ*� = 1

*μr*� ,

, *Fr*�(*r*) =

satisfy the cloaking coordinate transformation[69, 91], and they are *μr*� = *r*�

*γr*� = *γθ*� = *ω*0/100. *ωpz*� = *ωpr*� = *ωpθ*� = 4*ω*0, *Fz*�(*r*)=(*εz*� − 1)

, *Fθ*�(*r*)=(*μθ*� − 1)

where the light rays are forced to "diffract" around the cloaked area.

processes which we will discuss in another section [92].

*εz*� = ( *R*2 *R*2−*R*1

(*μr*� − 1)

before.

the related topics.

**7. Summary**

)2 *r*� −*R*1

> 0 )2+*ω*2 0*γ*2 *r*�

(*ω*2 *ar*�−*ω*2

(*ω*2 *ar*�−*ω*2 0 )*ω*2 *pr*�

[44] M. Bass, Handbook of Optics, Vol. II (McGraw-Hill, Inc., New York, 1995); W. Su, J. Li

The Group Velocity Picture of Metamaterial Systems 89

[46] E. Cubukcu, K. Aydin, E. Ozbay, S. Foteinopoulou, and C. M. Soukoulis Phys. Rev.Lett. 91, 207401 (2003); E. Cubukcu, K. Aydin, E. Ozbay, S. Foteinopoulou, and C.M.

[47] J. B. Pendry, Phys. Rev. Lett. 91, 099701 (2003);D. R. Smith, D. Schurig, M. Rosenbluth, S. Schultz, S. A. Ramakrishna, and J. B. Pendry, Appl. Phys. Lett. 82, 1506 (2003).

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

[50] X. S. Rao and C. K. Ong, Phys. Rev. B 68, 113103 (2003); X. S. Rao and C. K. Ong, Phys.

[52] Lei Zhou, C. T. Chan, Appl. Phys. Lett. 86, 101104 (2005), Y. Zhang, T. M. Grzegorczyk

[56] P. F. Loschialpo, D. L. Smith, D. W. Forester, F. J. Rachford, and J. Schelleng, Phys. Rev.

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

[60] Ilya V. Shadrivov, Andrey A. Sukhorukov, and Yuri S. Kivshar, Phys. Rev. E 67, 057602

[61] B. E. A. Saleh, M. C. Teich, *Fundamentals of Photonics* (John Wiley & Sons, New York,

[62] L. Mandel and E. Wolf *Optical Coherence and Quantum Optics* (Cambridge University, Cambridge, 1995); Marlan O. Scully and M. Suhail Zubairy *Quantum Optics*,

[64] In our source frequency range, the index range is about −1 − *i*0.0029 ± (0.006 + *i*10−6),

[65] The "group velocity" is not a well-defined value if the working frequency *ω*0 is near the

[69] S. A. Cummer, B.-I Popa, D. Schurig, and David R. Smith, Phys. Rev. E 74, 036621 (2006). [70] W. Cai, U. K. Chettiar, A. V. Kildishev, and V. M. Shalaev, Nat. Photonics 1, 224 (2007). [71] W. Cai, U. K. Chettiar, A. V. Kildishev, V. M. Shalaev, and G. W. Milton, Appl. Phys. Lett.

[72] D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R.

[73] H. Chen, B.-I. Wu, B. Zhang, and J. A. Kong, Phys. Rev. Lett. 99, 063903 (2007).

[59] Nader Engheta, IEEE Antennas and Wireless Propagation Lett. 1, 10, 2002

(2003); A. C. Peacock and N. G. R. Broderick, 11, 2502 (2003).

(2003); J. B. Pendry and D. R. Smith, Phys. Rev. Lett. 90, 029703 (2003).

and N. Xu, J. Biotechnol. 105, 165 (2003). [45] T. Tajima, Eur. Phys. J. D 55, 519 (2009).

Soukoulis, Nature (London) 423, 604 (2003).

[48] G. Gomez-Santos, Phys. Rev. Lett. 90, 077401 (2003).

[51] Michael W. Feise, Yuri S. Kivshar, Phys. Lett. A 334 326 (2005)

[53] L. Chen, S. He and L. Shen, Phys. Rev. Lett. 92, 107404(2004). [54] R. W. Ziolkowski and E. Heyman, Phys. Rev. E 64, 056625 (2001).

Rev. E 68, 067601 (2003).

E 67, 025602(R) (2003).

(2003)

1991)

and J. A. Kong, PIER 35,271, (2002)

[55] S. A. Cummer, Appl. Phys. Lett. 82, 1503 (2003).

[57] R. Merlin, Appl. Phys. Lett. 84, 1290 (2004).

(Cambridge University, Cambridge, 1997)

[66] Xunya Jiang *et. al.* unpublished. [67] U. Leonhardt, Science 312, 1777 (2006).

Smith, Science 314, 977 (2006).

91, 111105 (2007).

[63] X. Jiang and C. M. Soukoulis, Phys. Rev. Lett. 85, 70 (2000)

so the focal length defference and reflection are very small.

[68] J. B. Pendry, D. Schurig, and D. R. Smith, Science 312, 1780 (2006).

resonant frequency *ωa* of the NIM. But the GRT is still well-defined.

32 Will-be-set-by-IN-TECH

[15] J. B. Pendry, Opt. Express 11, 755 (2003); Zhaowei Liu, Hyesog Lee, Yi Xiong, Cheng Sun, and Xiang Zhang, Science 315, 1686 (2007); Igor I. Smolyaninov, Yu-Ju Hung, and

[16] Z. Jacob, L. V. Alekseyev, E. Narimanov, Opt. Express 14,8247 (2006); A. Salandrino, N. Engheta, Phys. Rev. B 74, 075103 (2006); Zubin Jacob, Leonid V. Alekseyev, and Evgenii Narimanov, J. Opt. Soc. Am. A 24,A52 (2007); J. B. Pendry, S. A. Ramakrishna, J. Phys.

[17] Xin Li,Zixian Liang,Xiaohan Liu,Xunya Jiang, and Jian Zi, Appl. Phys. Lett. 93, 171111

[19] P. W. Milonni, *Fast Light, Slow Light and Left-Handed Light* Institute of Physics, Great

[22] V. A. Podolskiy, and E. E. Narimanov, Phys. Rev. B 71, 201101 (2005); R. Wangberg, J. Elser, E. E. Narimanov, and V. A. Podolskiy, J. Opt. Soc. Am. B. 23, 498 (2006); A. A.

[23] A. Taflove, *Computational Electrodynamics: The Finite-Difference Time- Domain Method*

[26] Jie Yao,Zhaowei Liu,Yongmin Liu, Yuan Wang,Cheng Sun,Guy Bartal, Angelica M.

[31] N. Ashcroft, N. Mermin, Solid State Physics (Harcourt College Publishers, New York, 1976) ; H. Bruus and K. Flensberg, Many-body quantum theory in condensed matter

[33] R. Feynman, QED: the strange theory of light and matter (Princeton university press,

[39] R. G. Hunsperger, Integrated Optics: Theory and Technology (Springer-Verlag, New

[42] J. Kong, Electromagnetic Wave Theory, Electromagnetic Wave Theory

[20] As for the explanation of the static solution of the HI it'll be given in another paper.

[18] Hailu Luo, Wei Hu and Zhongzhou Ren, Europhys. Lett. 74 1081(2006).

[21] David R. Smith and Norman Kroll, Phys. Rev. Lett. 85, 2933(2000).

Govyadinov, V. A. Podolskiy, Phys. Rev. B 73,115108(2006).

[25] N. Fang, H. Lee, C. Sun, and X. Zhang, Science 308, 534 (2005).

Stacy and Xiang Zhang, Science 321, 930 (2008).

physics (Oxford university press, Oxford, 2004).

[35] Y. Ding, A. Kaplan, Phys. Rev. Lett. 63, 2725 (1989). [36] B. King, A. Piazza and C. Keitel, Nat. photon. 4, 92 (2010).

[37] H. Gies, Phys. Rev. D 61, 085021 (2000).

(Wiley-Interscience, New York, 1990).

[24] Alessandro Salandrino and Nader Engheta, Phys. Rev. B 74,115108(2006).

[27] J. B. Pendry and A. J. Holden and W. J. Stewart, Phys. Rev. Lett. 76, 4773(1996). [28] A. N. Grigorenko, A. K. Geim, H. F. Gleeson,*et al*, Nature, 438 335(2005)

[29] G. Dolling, M. Wegener, C. M. Soukoulis, and S. Linden, Opt. Lett. 32, 53 (2007). [30] D. Perkins, Introduction to high energy physics, 4th edition (Cambridge university

[32] J. Goodman, Statistical Optics (John Wiley & Sons, Inc., New York, 1985).

[34] R. Ruffini, G. Vereshchagin and S. Xue, Phys. Rep. 487, 1 (2010).

[38] J. S Heyl and L. Hernquist, J. Phys. A: Math. Gen. 30, 6485 (1997).

[40] B.Barwick, D.Flannigan, and A.H. Zewail, Nature 462, 902 (2009). [41] C. Carniglia and L. Mandel, J. Opt. Soc. Am. 61, 1035 (1971).

[43] J. Lawall and E. Kessler , Rev. Sci. Instrum. 71, 2669 (2000).

Christopher C. Davis, Science 315, 1699 (2007).

Condens. Matter 14,8463(2002).

Britain, CRC Press, 2004

(Artech House, Boston, 1995).

press, Cambridge, 2000).

Princeton, 1985).

York, 1985).

(2008).

**1. Introduction**

**4**

*USA*

Elena Semouchkina

*Michigan Technological University*

**Formation of Coherent Multi-Element**

Employment of metamaterials in cloaking devices, which could make concealed objects invisible, is based on the capability of resonators used as metamaterial "atoms" to provide effective material parameters, ranging from any positive to any negative values, in the vicinity of the resonance frequency. Metamaterials comprising layers with variously sized resonators (Schurig et al., 2006) or layers with variable density of identical resonators (Semouchkina et al., 2010) have been proposed to obtain the desired spatial distribution of the medium parameters, for example, radial variation of the effective permittivity or permeability in the shell prescribed by the transformation optics relations for cloaking cylindrically shaped objects (Leonhardt, 2006; Pendry et al., 2006). Such approach could be justified, if resonators contribute to the value of the effective parameters by the same way as polarized atoms contribute to the value of the polarization of a dielectric medium, i.e., if the responses of individual atoms are additive so that the medium response exhibits the same dispersion as the response of an individual atom. Then the material can be described by the effective medium theory (EMT), which was adopted for metamaterails description at the onset of metamaterial studies (Pendry et al., 1999; Pendry & Smith, 2003; Smith & Pendry, 2004; Smith et al., 2006). According to the EMT, a polar dielectric material could be represented by a set of parallel dipoles switching their directions by 1800 in dependence on the phase of the external field. Resonance modes in the "atoms" of metamaterials could also be represented by equivalent dipoles and, similarly to atoms in dielectrics, these dipoles are expected to respond accordingly, i.e. to follow the phase of the external field and demonstrate a coherent response within the half wavelength of incident radiation. However, observation of the performance of the infrared invisibility cloak designed from glass resonators (Semouchkina et al., 2010) has shown that obtaining a coherent response of the cloak structure presents a serious problem, if resonators in the metamaterail are electromagnetically coupled. The problem of coupling between resonators in metamaterials,

It is well known that metamaterials, as a rule, are designed as periodic structures of closepacked resonators. When the dimensions of resonators are ten times smaller than the wavelength, the EMT is considered to be applicable, so that the waves "do not see" the atoms and propagate in the metamaterial as in a homogenized, i.e. in a uniform medium. Such homogenization is supposed to make metamaterials essentially different from

however, remains to be largely ignored in the literature.

**Resonance States in Metamaterials**

### **Formation of Coherent Multi-Element Resonance States in Metamaterials**

Elena Semouchkina *Michigan Technological University USA*

#### **1. Introduction**

34 Will-be-set-by-IN-TECH

90 Metamaterial

[79] A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference*

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

[81] H. Chen, Z. Liang, P. Yao, X. Jiang, H. Ma, and C. T. Chan, Phys. Rev. B 76, 241104 (2007). [82] R. L. Fante, and M. T. McCormack, IEEE Trans. Antennas Propag. 30, 1443(1968). [83] M. Kerker, J. Opt. Soc. Am. 65, 376(1975). and D. R. Nature Photonics 1, 224 (2007).

[84] Hongsheng Chen, Bae-Ian Wu, Baile Zhang, and Jin Au Kong, Phys. Rev. Lett. 99, 063903

[85] Zhichao Ruan, Min Yan, Curtis W. Neff, and Min Qiu, Phys. Rev. Lett. 99, 113903 (2007).

[91] When *r*� → *R*1, *μθ*� will tend to infinite. In order that it can be realizable in our numerical

[94] A. Greenleaf, Y. Kurylev, M. Lassas, G. Uhlmann, http:

[95] L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media (Pergamon

[96] HY Chen and CT Chan, the brief report for the detail of the causality constraint. [97] T. Koschny, M. Kafesaki, E. N. Economou and C. M. Soukoulis, Phys. Rev. Lett. 93,

[74] Z. Ruan, M. Yan, C. W. Neff, and M. Qiu, Phys. Rev. Lett. 99, 113903 (2007).

[77] X. Jiang, W. Han, P. Yao, and W. Li, Appl. Phys. Lett. 89, 221102 (2006). [78] P. Yao, W. Li, S. Feng, and X. Jiang, Opt. Express 14, 12295 (2006).

*Time-Domain Method*, 2nd ed. (Artech House, Boston, 2000).

Graeme W. Milton, Appl. Phys. Lett. 91, 111105 (2007).

[86] Huanyang Chen, Xunya Jiang, C. T. Chan, arXiv:0707.1126v2.

[88] H. Chen and C. T. Chan, Appl. Phys. Lett. 90, 241105 (2007). [89] Ulf Leonhardt and Thomas G Philbin, New J. Phys. 8, 247(2006). [90] J. B. Pendry and D. R. Smith, Phys. Rev. Lett. 90, 029703 (2003).

simulation, we limit its maximum value to 103.

//www.arXiv:/0703059[math-ph] (2007).

Press, 1975), ch. 11, pp. 315-321

107402 (2004).

[87] D. Schurig, J. B. Pendry, and D. R. Smith, Optics Express 14, 9794(2006).

[92] Zixian Liang, Peijun Yao, Xunya Jiang, and Xiaowei Sun, unpublished. [93] G. W. Milton, M. Briane, and J. R. Willis, New J. Phys. 8, 248 (2006).

[75] U. Leonhardt, New J. Phys. 8, 118 (2006). [76] S. A. Cummer, Appl. Phys. Lett. 82, 2008 (2003).

unpublished.

(2007).

Employment of metamaterials in cloaking devices, which could make concealed objects invisible, is based on the capability of resonators used as metamaterial "atoms" to provide effective material parameters, ranging from any positive to any negative values, in the vicinity of the resonance frequency. Metamaterials comprising layers with variously sized resonators (Schurig et al., 2006) or layers with variable density of identical resonators (Semouchkina et al., 2010) have been proposed to obtain the desired spatial distribution of the medium parameters, for example, radial variation of the effective permittivity or permeability in the shell prescribed by the transformation optics relations for cloaking cylindrically shaped objects (Leonhardt, 2006; Pendry et al., 2006). Such approach could be justified, if resonators contribute to the value of the effective parameters by the same way as polarized atoms contribute to the value of the polarization of a dielectric medium, i.e., if the responses of individual atoms are additive so that the medium response exhibits the same dispersion as the response of an individual atom. Then the material can be described by the effective medium theory (EMT), which was adopted for metamaterails description at the onset of metamaterial studies (Pendry et al., 1999; Pendry & Smith, 2003; Smith & Pendry, 2004; Smith et al., 2006). According to the EMT, a polar dielectric material could be represented by a set of parallel dipoles switching their directions by 1800 in dependence on the phase of the external field. Resonance modes in the "atoms" of metamaterials could also be represented by equivalent dipoles and, similarly to atoms in dielectrics, these dipoles are expected to respond accordingly, i.e. to follow the phase of the external field and demonstrate a coherent response within the half wavelength of incident radiation. However, observation of the performance of the infrared invisibility cloak designed from glass resonators (Semouchkina et al., 2010) has shown that obtaining a coherent response of the cloak structure presents a serious problem, if resonators in the metamaterail are electromagnetically coupled. The problem of coupling between resonators in metamaterials, however, remains to be largely ignored in the literature.

It is well known that metamaterials, as a rule, are designed as periodic structures of closepacked resonators. When the dimensions of resonators are ten times smaller than the wavelength, the EMT is considered to be applicable, so that the waves "do not see" the atoms and propagate in the metamaterial as in a homogenized, i.e. in a uniform medium. Such homogenization is supposed to make metamaterials essentially different from

Formation of Coherent Multi-Element Resonance States in Metamaterials 93

Fig. 1. Cross-section of the metamaterial block depicting SRR sections (longer strips) and cut wires (shorter strips) placed normally to the figure plane (left figure); and distributions of magnetic field amplitude (upper row) and phase (lower row) in the same cross-section at four frequencies within the metamaterial transmission band. Bright spots in the upper row

This Chapter addresses the listed above underexplored metamaterial problems and describes the effects caused by coupling and splitting of resonances in two types of metamaterial multi-resonator structures, i.e. those composed of SRRs and of dielectric resonators. The results presented here demonstrate that EM responses provided by these structures cannot be reduced to the responses of their constituent single resonators as it is assumed at application of the EMT. The Chapter also shows that splitting of resonances in arrays of dielectric resonators leads to the formation of a bandgap dividing two transmission zones, one at lower, and another one at higher frequncies. While at the gap frequencies elementary resonances in arrays look sporadically flashing at various multiple frequencies, the lower frequency transmission zone is characterized by overlapping of the resonance fields and slow wave propagation. The transmission zone at higher frequencies can be characterized by a coherent response of elementary resonators within halfwavelength distances and by superluminal phase velocities of propagating waves. The obtained results bridge the gap between the properties of metamaterials and photonic

**2. Resonance splitting in metamaterial arrays composed of SRRs**

Investigation of the resonance splitting phenomena in SRR, which is presented in the following sub-sections, has been performed on fragments of SRR arrays quite similar to the arrays used in the design of the microwave invisibility cloak (Schurig et al., 2006). This cloaking shell was formed from concentric arrays of unit cells each including one SRR. The size of the sell was of 3.33 mm x 3.33 mm in the SRR plane and 3.18 mm normal to this plane with the SRR located in the central position. The SRRs were placed on the substrate with the

correspond to the areas with highest field magnitude.

crystals employing Mie resonances.

**2.1 Specific features of SRR characterization**

photonic crystals (Pendry & Smith, 2003). It is worth noting that the lattice parameter was not thought to be critical for the EMT application, so that the only ultimate requirement for metamaterials in this attitude is that the lattice parameter should not be equal to the half wavelength of incident radiation, which is typical for photonic crystals. Therefore, for example, the metamaterial for the microwave cloak operating at 8.5 GHz (Schurig et al., 2006), was composed of split ring resonators (SRR) with the planar dimensions of 3 mm x 3 mm packed in concentric arrays with the inter-resonators separation of 0.17 mm, which is negligible compared to the resonator size. Such close packing of SRRs causes questions about possible effects of interaction between resonators and elementary resonance splitting that could deteriorate the EMT applicability.

Our earlier studies (Semouchkina et al., 2004, 2005) indeed pointed out at strong splitting of elementary resonances in close-packed resonator arrays of conventional metamaterials consisting of metal elements. At the modeling of an extended finite metamaterial block, which consisted of 36 SRRs and 12 elongated metal strips arranged in unit cells quite similar to typical arrangements employed in conventional metamaterials, it was shown that at the frequencies corresponding to the metamaterial transmission band, elementary resonances were essentially coupled and integrated in 3D networks, an example of which is presented in Fig. 1. As a sequence of coupling and splitting, most resonators responded resonantly at multiple frequencies and formed specific accociations at different frequencies in the transmission band, making the metamaterial essentially nonuniform with the character of inhomogenities changing in dependence on frequency. No regular wave front propagation across the sample was observed, so that the energy transfer occurred rather due to the hopping mechanism similar to that thought to be inherent in the CROWs (coupled resonator optical waveguides (Yariv et al., 1999).

It is logical to suggest that phenomena similar to the mentioned above should take place in any close-packed metamaterial structures including those used in the microwave cloak described in (Schurig et al., 2006). The reason why these phenomena have not been revealed by the authors can be seen in the employment of simulation models, in which real multiresonator cloak structures were replaced by layered material structures with prescribed values of the effective permeability for each layer.

An additional reason, why coupling phenomena in metamaterials remain underexplored is seen in the typical approach to consider metamaterial properties to be identical to the properties of a single cell and to ignore interaction between cells. This approach is based on the results of unit cell simulations which could not reveal resonance splitting even at small distances between neighboring resonators defined by the dimensions of the unit cell, because of employment of plane wave excitation source at normal wave incidence and periodic boundary conditions (PBC) at the cell boundaries normal to the direction of wave propagation. Such arangement, however, could only model an infinite 2D array of unit cells and not a 3D metamaterial medium. Our results indicate that drastic changes in electromagnetic (EM) metamaterial response could be observed at stacking 2D arrays in the direction of wave propagation to form a 3D structure, because incident wave has a different phase, when encounters each 2D array. Significant changes also occur in finite 2D arrays and in arrays with some distortions of periodicity that is unavoidable in practice.

photonic crystals (Pendry & Smith, 2003). It is worth noting that the lattice parameter was not thought to be critical for the EMT application, so that the only ultimate requirement for metamaterials in this attitude is that the lattice parameter should not be equal to the half wavelength of incident radiation, which is typical for photonic crystals. Therefore, for example, the metamaterial for the microwave cloak operating at 8.5 GHz (Schurig et al., 2006), was composed of split ring resonators (SRR) with the planar dimensions of 3 mm x 3 mm packed in concentric arrays with the inter-resonators separation of 0.17 mm, which is negligible compared to the resonator size. Such close packing of SRRs causes questions about possible effects of interaction between resonators and elementary resonance splitting

Our earlier studies (Semouchkina et al., 2004, 2005) indeed pointed out at strong splitting of elementary resonances in close-packed resonator arrays of conventional metamaterials consisting of metal elements. At the modeling of an extended finite metamaterial block, which consisted of 36 SRRs and 12 elongated metal strips arranged in unit cells quite similar to typical arrangements employed in conventional metamaterials, it was shown that at the frequencies corresponding to the metamaterial transmission band, elementary resonances were essentially coupled and integrated in 3D networks, an example of which is presented in Fig. 1. As a sequence of coupling and splitting, most resonators responded resonantly at multiple frequencies and formed specific accociations at different frequencies in the transmission band, making the metamaterial essentially nonuniform with the character of inhomogenities changing in dependence on frequency. No regular wave front propagation across the sample was observed, so that the energy transfer occurred rather due to the hopping mechanism similar to that thought to be inherent in the CROWs (coupled resonator

It is logical to suggest that phenomena similar to the mentioned above should take place in any close-packed metamaterial structures including those used in the microwave cloak described in (Schurig et al., 2006). The reason why these phenomena have not been revealed by the authors can be seen in the employment of simulation models, in which real multiresonator cloak structures were replaced by layered material structures with prescribed

An additional reason, why coupling phenomena in metamaterials remain underexplored is seen in the typical approach to consider metamaterial properties to be identical to the properties of a single cell and to ignore interaction between cells. This approach is based on the results of unit cell simulations which could not reveal resonance splitting even at small distances between neighboring resonators defined by the dimensions of the unit cell, because of employment of plane wave excitation source at normal wave incidence and periodic boundary conditions (PBC) at the cell boundaries normal to the direction of wave propagation. Such arangement, however, could only model an infinite 2D array of unit cells and not a 3D metamaterial medium. Our results indicate that drastic changes in electromagnetic (EM) metamaterial response could be observed at stacking 2D arrays in the direction of wave propagation to form a 3D structure, because incident wave has a different phase, when encounters each 2D array. Significant changes also occur in finite 2D arrays and in arrays with some distortions of periodicity that is unavoidable in

that could deteriorate the EMT applicability.

optical waveguides (Yariv et al., 1999).

practice.

values of the effective permeability for each layer.

Fig. 1. Cross-section of the metamaterial block depicting SRR sections (longer strips) and cut wires (shorter strips) placed normally to the figure plane (left figure); and distributions of magnetic field amplitude (upper row) and phase (lower row) in the same cross-section at four frequencies within the metamaterial transmission band. Bright spots in the upper row correspond to the areas with highest field magnitude.

This Chapter addresses the listed above underexplored metamaterial problems and describes the effects caused by coupling and splitting of resonances in two types of metamaterial multi-resonator structures, i.e. those composed of SRRs and of dielectric resonators. The results presented here demonstrate that EM responses provided by these structures cannot be reduced to the responses of their constituent single resonators as it is assumed at application of the EMT. The Chapter also shows that splitting of resonances in arrays of dielectric resonators leads to the formation of a bandgap dividing two transmission zones, one at lower, and another one at higher frequncies. While at the gap frequencies elementary resonances in arrays look sporadically flashing at various multiple frequencies, the lower frequency transmission zone is characterized by overlapping of the resonance fields and slow wave propagation. The transmission zone at higher frequencies can be characterized by a coherent response of elementary resonators within halfwavelength distances and by superluminal phase velocities of propagating waves. The obtained results bridge the gap between the properties of metamaterials and photonic crystals employing Mie resonances.

### **2. Resonance splitting in metamaterial arrays composed of SRRs**

#### **2.1 Specific features of SRR characterization**

Investigation of the resonance splitting phenomena in SRR, which is presented in the following sub-sections, has been performed on fragments of SRR arrays quite similar to the arrays used in the design of the microwave invisibility cloak (Schurig et al., 2006). This cloaking shell was formed from concentric arrays of unit cells each including one SRR. The size of the sell was of 3.33 mm x 3.33 mm in the SRR plane and 3.18 mm normal to this plane with the SRR located in the central position. The SRRs were placed on the substrate with the

Formation of Coherent Multi-Element Resonance States in Metamaterials 95

been characterized. The values of *s* and *r* were taken very close to those characteristic for the SRRs of the 1st, 2nd, and 3rd concentric arrays of the microwave cloak described in (Schurig et al., 2006). Small corrections, however, have been introduced in order to provide fitting of the chosen *s* and *r* to the least square dependences which were built to express interrelation between these parameters of SRRs and the effective permeability found in (Schurig et al.,

mm mm GHz GHz GHz GHz

7.761

9.125

9.061

As seen from Table 1, numerical experiments with plane wave incidence and with the standard waveguide gave quite close results for *res f* , while in experiments with unit cells used in (Schurig et al., 2006) this frequency appeared to be essentially lower both for SRRs placed in air and on a substrate. The resulting S-parameter spectra in the latter simulations appeared distorted in comparison to the spectra obtained in the standard waveguide. Similar distortions of resonances are known to be observed at too close placement of metal objects or the domain boundaries to the resonator. Therefore, the simulations have been repeated for computational volumes of increased dimensions. Gradual increase of the volume caused better and better fitting of the simulation results to the data obtained in two other sets of numerical experiments. Saturation of changes and coincidence of the data was observed at the distances between resonator and domain boundaries exceeding 1.833 mm in the SRR plane and 4.183 mm in the normal direction. Apparently these dimensions characterize the "halo" of the resonance fields around the SRRs, illustrated by the field

It is obvious that the space occupied by the SRR field "halo" substantially exceeds the volume of unit cells used in (Schurig et al., 2006). As seen in Figs. 2a and 2c, the resonance magnetic and electric fields in the "halo" are not vanishing at the distances from the SRR

*free space fres*

7.673 7.657

7.745 7.19 7.76

*unit cell f res*

7.125

9.16 9.119 8.53 9.17

7.73 7.67 7.73

9.075

*waveguide f res*

Transient solver Frequency domain

8.448

9.03 9.068

7.027 7.673

8.49 9.185

solver

*halo cell f res*

2006) for unit cells with SRRs of proper dimensions.

*r*

*substrate*

yes

no

yes

no

no

Table 1. Parameters and resonance frequencies of various SRRs

yes

1 1.654 0.26

2 1.677 0.254

patterns presented in Figs. 2a and 2c.

1.716 0.244

*s*

N

Table 1

3

thickness of 381 microns and the relative permittivity of 2.33. In the simulations described in the next sub-section, the SRRs with identical parameters (Fig. 2a) were used to build arrays of unit cells similar to various fragments of the cloak described in (Schurig et al., 2006). The only difference between the investigated metamaterial fragments and those used in the publication was that fragments in our case were flat. However, as the radiuses of concentric arrays in (Schurig et al., 2006) were relatively big compared to the cell size, flat arrays with the lengths of up to six cells were a good approaximation of fractions of concentric arrays in the cloak metamaterial. The propagation vector (k-vector) of the incident wave in our studies was always directed along the array rows in the plane of the SRR location, while magnetic field was directed normally to the SRR planes. Therefore, the investigated arrays had to perform similarly to the cloak fragments located close to the diameter of the shell normal to the k-vector of the incident wave. First, the SRR responses in air and at the placement of the resonators on the dielectric substrate have been investigated. Numerical simulations were performed for three arrangements: 1) at the plane wave incidence on the object, 2) at the placement of the SRR in the waveguide of standard dimensions with perfect electric conductor (PEC) walls, and 3) at application of waveguide ports to the standard unit cell and/or to the unit cell increased to incorporate the resonance field "halos" (the following section describes how the "halo" dimensions were found).

Fig. 2. Resonance in a single SRR placed in the waveguide WR137: (a) measured and simulated transmission spectra S21; (b) arrangement for the measurements; (c) magnetic and (d) electric field patterns.

All simulations were conducted by using the commercial full-wave software package CST Microwave Studio. For the two first arrangements the Transient Solver has been used, while for the third one – the Frequency Domain Solver. Numerical experiments with waveguide ports allowed for obtaining the scattering parameter spectra convenient for controlling the resonance frequency *res f* , while in experiments with plane wave incidence, the resonance was detected by monitoring the signal of H-field probes placed at the expected maximum field locations in the resonators. Three types of SRRs that differed by the lengths of their inner slot forming strips *s* and the radii of the corners *r,* as summarized in Table 1, have

thickness of 381 microns and the relative permittivity of 2.33. In the simulations described in the next sub-section, the SRRs with identical parameters (Fig. 2a) were used to build arrays of unit cells similar to various fragments of the cloak described in (Schurig et al., 2006). The only difference between the investigated metamaterial fragments and those used in the publication was that fragments in our case were flat. However, as the radiuses of concentric arrays in (Schurig et al., 2006) were relatively big compared to the cell size, flat arrays with the lengths of up to six cells were a good approaximation of fractions of concentric arrays in the cloak metamaterial. The propagation vector (k-vector) of the incident wave in our studies was always directed along the array rows in the plane of the SRR location, while magnetic field was directed normally to the SRR planes. Therefore, the investigated arrays had to perform similarly to the cloak fragments located close to the diameter of the shell normal to the k-vector of the incident wave. First, the SRR responses in air and at the placement of the resonators on the dielectric substrate have been investigated. Numerical simulations were performed for three arrangements: 1) at the plane wave incidence on the object, 2) at the placement of the SRR in the waveguide of standard dimensions with perfect electric conductor (PEC) walls, and 3) at application of waveguide ports to the standard unit cell and/or to the unit cell increased to incorporate the resonance field "halos" (the

k H

E

(b)

(b)

simulated measured

Fig. 2. Resonance in a single SRR placed in the waveguide WR137: (a) measured and simulated transmission spectra S21; (b) arrangement for the measurements; (c) magnetic and

All simulations were conducted by using the commercial full-wave software package CST Microwave Studio. For the two first arrangements the Transient Solver has been used, while for the third one – the Frequency Domain Solver. Numerical experiments with waveguide ports allowed for obtaining the scattering parameter spectra convenient for controlling the resonance frequency *res f* , while in experiments with plane wave incidence, the resonance was detected by monitoring the signal of H-field probes placed at the expected maximum field locations in the resonators. Three types of SRRs that differed by the lengths of their inner slot forming strips *s* and the radii of the corners *r,* as summarized in Table 1, have

following section describes how the "halo" dimensions were found).

6.5 7 7.5 8 8.5 9

Frequency (GHz)

(c) (d)

S21 (dB)

(d) electric field patterns.

(d) (a)

been characterized. The values of *s* and *r* were taken very close to those characteristic for the SRRs of the 1st, 2nd, and 3rd concentric arrays of the microwave cloak described in (Schurig et al., 2006). Small corrections, however, have been introduced in order to provide fitting of the chosen *s* and *r* to the least square dependences which were built to express interrelation between these parameters of SRRs and the effective permeability found in (Schurig et al., 2006) for unit cells with SRRs of proper dimensions.

Table 1. Parameters and resonance frequencies of various SRRs

As seen from Table 1, numerical experiments with plane wave incidence and with the standard waveguide gave quite close results for *res f* , while in experiments with unit cells used in (Schurig et al., 2006) this frequency appeared to be essentially lower both for SRRs placed in air and on a substrate. The resulting S-parameter spectra in the latter simulations appeared distorted in comparison to the spectra obtained in the standard waveguide. Similar distortions of resonances are known to be observed at too close placement of metal objects or the domain boundaries to the resonator. Therefore, the simulations have been repeated for computational volumes of increased dimensions. Gradual increase of the volume caused better and better fitting of the simulation results to the data obtained in two other sets of numerical experiments. Saturation of changes and coincidence of the data was observed at the distances between resonator and domain boundaries exceeding 1.833 mm in the SRR plane and 4.183 mm in the normal direction. Apparently these dimensions characterize the "halo" of the resonance fields around the SRRs, illustrated by the field patterns presented in Figs. 2a and 2c.

It is obvious that the space occupied by the SRR field "halo" substantially exceeds the volume of unit cells used in (Schurig et al., 2006). As seen in Figs. 2a and 2c, the resonance magnetic and electric fields in the "halo" are not vanishing at the distances from the SRR

Formation of Coherent Multi-Element Resonance States in Metamaterials 97

7.1 7.3 7.5 7.7 7.9 8.1 8.3 8.5

2-SRR 3-SRR 4-SRR 5-SRR

**degrees**

**degrees**

Frequency (GHz)

1940, 2020, 15.60, 25.10, 25.20

Fig. 3. Resonances in linear arrays arranged in wave propagation direction in the waveguide WR137: (a) S21; (b) coherent response of a 5-SRR array at 7.73 GHz; and (c) assymmetric

Much stronger splitting of SRR resonances has been revealed in the SRR columns (Fig. 4) similar to those used in concentric arrays of the SRRs in (Schurig et al., 2006). Placement of such columns in the waveguide WR137 caused some concerns, since at about 1 cm height of the columns their resonance fields could interact with the waveguide top and bottom walls that could cause deterioration of the resonances. Therefore, the simulations for these columns were initially performed in free space at plane wave incidence. As such simulations did not provide an opportunity for deriving the transmission spectra, the resonance field in the columns was monitored by using H-field probes placed in geometrically identical points of the SRR cross-sections (Fig. 4a). As seen from Fig.4e the probes revealed three resonances covering more than 1 GHz of the frequency range (i.e. almost five times wider resonance band than that observed for the linear arrays of SRRs arranged along the direction of the incident wave k—vector). It is worth mentioning that each resonator in the column responded resonantly at all split frequencies while the integrated modes at these frequencies were quite different. As seen from Fig. 4b, at the lower frequency resonance sin-phase oscillations were observed in the two upper resonators, while the third resonator oscillated with a smaller magnitude and with opposite phase. At the higher frequency resonance, the two lower resonators demonstrated sin-phase oscillations while the third resonator oscillated with opposite phase (Fig. 4d). Resonance at the median frequency was supported mainly by sin-phase oscillations in the upper and the lower resonators of the column (Fig. 4c). The third resonator located in the middle of the column did not show strong

1130, 1110, 1110, 1100, 1110

S21 (dB)

(b)

(c)

response at 7.95 GHz.

0

(a)

corresponding to the boundaries of the used unit cells, therefore, the resonance formation and the resonance frequencies *res f* were affected in respective simulations of such unit cells. In comparison, field magnitudes at the boundaries of unit cells with the size equal to the "halo" size were found to be ten times smaller than those of the cells used in (Schurig et al., 2006) that explains the consistency of the data obtained for the halo-sized cells and the results of other numerical experiments. Only at building the shell of such cells, i.e. at the placement of resonators in arrays of the metamaterial at distances, when "halos" are not overlapped, interaction between resonators could be omitted from consideration. In contrast, the distance between resonators in concentric arrays in (Schurig et al., 2006) did not exceed 0.333 mm, while the distance between arrays (i.e. along the normal to the plane of resonators) was only 3.183 mm.

The data presented in Table 1 do not confirm the expectations for *res f* to be of about 8.5 GHz as reported in (Schurig et al., 2006) for both the resonance frequency of basic SRRs and the frequency of the cloak response. Resonances at frequencies close to 8.5 GHz were obtained in our experiments only for the unit cells used in (Schurig et al., 2006) at the placement of SRRs in air. It might be that the dielectric substrate was omitted in models used for simulations in (Schurig et al., 2006). The obtained discrepancy does not seem to be critical for the studies of the splitting phenomena in SRR-based metamaterials employed in (Schurig et al., 2006), although it should be taken into account at verification of the cloaking effect, since this effect should not be expected at 8.5 GHz, if the metamaterial perform in accordance to the EMT predictions.

#### **2.2 Resonances in SRR arrays representing fragments of the microwave cloak**

In order to additionally verify the results of simulations, individual SRRs and various multielement arrays of SRRs with the same dimensions as in (Schurig et al., 2006) have been fabricated, and their EM response has been measured in the waveguide WR137 operating at the single TE10 mode in the range (5-10) GHz. In experiments, the waveguide WR137 loaded with SRR arrays was connected to the PNA-L Network Analyzer N5230A to measure the transmission spectrum. In order to avoid EM interaction between the arrays and waveguide walls, resonators or arrays were supported by a styrofoam layer inserted into the center region of the waveguide. As seen from Fig. 2d, the resonance in a single resonator placed in the waveguide was observed at f=7.78 GHz and not at 8.5 GHz, i.e. quite close to the results of our simulations for the resonator with the same dimensions.

Next, linear arrays of SRRs arranged along the direction of k-vector of the incident wave have been investigated. Fig. 3a shows that the arrays of 3-5 SRRs demonstrate strong splitting of the resonance, which increases with the number of resonators in the array. It is interesting to note, that instead of multiple dips in the transmission spectrum correlated with the quantity of resonators in the array, only two well expressed drops in the transmission spectra have been observed. The first drop is accompanied by the coherent response of resonators (Fig. 3b), while at the higher frequency drop the resonators in the line respond with opposite phases in two halves of the array. The phase patterns of the wave propagating in free space inside the waveguide at both resonances show strong disturbances of the wave transmission through the waveguide (the curve between bright and dark parts of the image divides the areas corresponding to the half wavelengths of the propagating waves).

corresponding to the boundaries of the used unit cells, therefore, the resonance formation and the resonance frequencies *res f* were affected in respective simulations of such unit cells. In comparison, field magnitudes at the boundaries of unit cells with the size equal to the "halo" size were found to be ten times smaller than those of the cells used in (Schurig et al., 2006) that explains the consistency of the data obtained for the halo-sized cells and the results of other numerical experiments. Only at building the shell of such cells, i.e. at the placement of resonators in arrays of the metamaterial at distances, when "halos" are not overlapped, interaction between resonators could be omitted from consideration. In contrast, the distance between resonators in concentric arrays in (Schurig et al., 2006) did not exceed 0.333 mm, while the distance between arrays (i.e. along the normal to the plane of

The data presented in Table 1 do not confirm the expectations for *res f* to be of about 8.5 GHz as reported in (Schurig et al., 2006) for both the resonance frequency of basic SRRs and the frequency of the cloak response. Resonances at frequencies close to 8.5 GHz were obtained in our experiments only for the unit cells used in (Schurig et al., 2006) at the placement of SRRs in air. It might be that the dielectric substrate was omitted in models used for simulations in (Schurig et al., 2006). The obtained discrepancy does not seem to be critical for the studies of the splitting phenomena in SRR-based metamaterials employed in (Schurig et al., 2006), although it should be taken into account at verification of the cloaking effect, since this effect should not be expected at 8.5 GHz, if the metamaterial perform in

**2.2 Resonances in SRR arrays representing fragments of the microwave cloak**

of our simulations for the resonator with the same dimensions.

corresponding to the half wavelengths of the propagating waves).

In order to additionally verify the results of simulations, individual SRRs and various multielement arrays of SRRs with the same dimensions as in (Schurig et al., 2006) have been fabricated, and their EM response has been measured in the waveguide WR137 operating at the single TE10 mode in the range (5-10) GHz. In experiments, the waveguide WR137 loaded with SRR arrays was connected to the PNA-L Network Analyzer N5230A to measure the transmission spectrum. In order to avoid EM interaction between the arrays and waveguide walls, resonators or arrays were supported by a styrofoam layer inserted into the center region of the waveguide. As seen from Fig. 2d, the resonance in a single resonator placed in the waveguide was observed at f=7.78 GHz and not at 8.5 GHz, i.e. quite close to the results

Next, linear arrays of SRRs arranged along the direction of k-vector of the incident wave have been investigated. Fig. 3a shows that the arrays of 3-5 SRRs demonstrate strong splitting of the resonance, which increases with the number of resonators in the array. It is interesting to note, that instead of multiple dips in the transmission spectrum correlated with the quantity of resonators in the array, only two well expressed drops in the transmission spectra have been observed. The first drop is accompanied by the coherent response of resonators (Fig. 3b), while at the higher frequency drop the resonators in the line respond with opposite phases in two halves of the array. The phase patterns of the wave propagating in free space inside the waveguide at both resonances show strong disturbances of the wave transmission through the waveguide (the curve between bright and dark parts of the image divides the areas

resonators) was only 3.183 mm.

accordance to the EMT predictions.

Fig. 3. Resonances in linear arrays arranged in wave propagation direction in the waveguide WR137: (a) S21; (b) coherent response of a 5-SRR array at 7.73 GHz; and (c) assymmetric response at 7.95 GHz.

Much stronger splitting of SRR resonances has been revealed in the SRR columns (Fig. 4) similar to those used in concentric arrays of the SRRs in (Schurig et al., 2006). Placement of such columns in the waveguide WR137 caused some concerns, since at about 1 cm height of the columns their resonance fields could interact with the waveguide top and bottom walls that could cause deterioration of the resonances. Therefore, the simulations for these columns were initially performed in free space at plane wave incidence. As such simulations did not provide an opportunity for deriving the transmission spectra, the resonance field in the columns was monitored by using H-field probes placed in geometrically identical points of the SRR cross-sections (Fig. 4a). As seen from Fig.4e the probes revealed three resonances covering more than 1 GHz of the frequency range (i.e. almost five times wider resonance band than that observed for the linear arrays of SRRs arranged along the direction of the incident wave k—vector). It is worth mentioning that each resonator in the column responded resonantly at all split frequencies while the integrated modes at these frequencies were quite different. As seen from Fig. 4b, at the lower frequency resonance sin-phase oscillations were observed in the two upper resonators, while the third resonator oscillated with a smaller magnitude and with opposite phase. At the higher frequency resonance, the two lower resonators demonstrated sin-phase oscillations while the third resonator oscillated with opposite phase (Fig. 4d). Resonance at the median frequency was supported mainly by sin-phase oscillations in the upper and the lower resonators of the column (Fig. 4c). The third resonator located in the middle of the column did not show strong

Formation of Coherent Multi-Element Resonance States in Metamaterials 99

responses with respect to the central row i.e. the resonators located in upper and lower rows responded coherently, even though neighbouring columns demonstrated opposite phases of oscillations. So, it could be said that every other column responded coherently (Fig. 6c).

90.40

90.40

1340 1340

> 4 5

13 14

15

7 8

9

6

6.5 7 7.5 8 8.5 9

6.5 7 7.5 8 8.5 9

column) for SRR pairs arranged: (a) "face-to-face" and (b) "back-to-back ".

6.95 GHz 7.7 GHz

frequency, GHz

Fig. 5. Geometry, S21 spectra, and patterns of H-field magnitude (3d column) and phase (4th

6.5 7.0 7.5 8.0 8.5 8.8 frequency, GHz

(b) (c) (d)8.15 GHz

Fig. 6. (a) The spectra of the probe signals; inset shows the SRR array and the H-field probe locations; (b)-(d) typical patterns of the resonance oscillations in the array sampled at 6.64,

Despite of the above mentioned concerns about possible interaction between multi-element arrays and waveguide walls at the measurements we have performed comparison of simulation and measurement results for the transmission spectra at the placement of the 3 x 5 arrays of SRRs in the waveguide WR137. A good matching of the results presented in Fig. 7 is obvious, although both S21 spectra demostrate a smaller quantity of the transmission

7.4 and 8.15 GHz (colour intensity shows the magnitude of resonance oscillations).

1 2 10 11

12

3

S21(dB)

*E H k*

(a)

(b)

0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0.0

(a)

Probe Magnitude in A/m

0

0

S21(dB)

oscillations, however, the phase of these oscillations was shifted by 1800 with respect to the oscillations of two other SRRs.

Fig. 4. (a) The model of the SRR column with H-field probes; (b)-(d) H-field resonance responses of the column at signal peak frequencies (colour intensity shows the magnitude of resonance; and (e) the spectra of probe signals.

It is interesting to note that shorted columns of 2 SRRs placed in either "face-to-face" or "back-to-back" arrangement did not demonstrate any splitting of the resonances. As seen from Fig. 5, the former arrangement provided for the unsplit resonance at almost the same frequency as the lower frequency resonance in the column of 3 SRRs, while the latter arrangement provided for the resonance at the frequency close to the frequency of upper resonance in the column of 3 SRRs. It means that the listed above resonances in a 3-SRR column are apparently related to integrated coherent resonances in couples of oppositely arranged SRRs. The coherentness is evident from the phase patterns presented in Fig. 5 on the right.

Essentially more complicated response with multiple resonances of various Q factors was demonstrated by arrays composed of several columns as, for example, the presented by the insert in Fig 6a array of 5 columns. Despite of the complexity, however, it was possible to distinguish three main groups of resonances comparable to resonances revealed in one SRR column (Fig. 4). In fact, the lower frequency group in the 3 x 5 SRR array demontsrated coherent resonances in upper rows of the array (Fig. 6b), while the higher frequency group demontsrated coherent resonances in the lower rows of the array (Fig. 6d). The group of resonances at the median frequencies could be characterized by a symmetry of DR

oscillations, however, the phase of these oscillations was shifted by 1800 with respect to the

(a) (b) (c) (d)

3

2

1

1 - 2 - 3 - 6.9 GHz 7.5 GHz 7.9 GHz

6.2 6.5 7.0 7.5 8.0 8.5 Frequency, GHz

Fig. 4. (a) The model of the SRR column with H-field probes; (b)-(d) H-field resonance responses of the column at signal peak frequencies (colour intensity shows the magnitude of

It is interesting to note that shorted columns of 2 SRRs placed in either "face-to-face" or "back-to-back" arrangement did not demonstrate any splitting of the resonances. As seen from Fig. 5, the former arrangement provided for the unsplit resonance at almost the same frequency as the lower frequency resonance in the column of 3 SRRs, while the latter arrangement provided for the resonance at the frequency close to the frequency of upper resonance in the column of 3 SRRs. It means that the listed above resonances in a 3-SRR column are apparently related to integrated coherent resonances in couples of oppositely arranged SRRs. The coherentness is evident from the phase patterns presented in Fig. 5 on

Essentially more complicated response with multiple resonances of various Q factors was demonstrated by arrays composed of several columns as, for example, the presented by the insert in Fig 6a array of 5 columns. Despite of the complexity, however, it was possible to distinguish three main groups of resonances comparable to resonances revealed in one SRR column (Fig. 4). In fact, the lower frequency group in the 3 x 5 SRR array demontsrated coherent resonances in upper rows of the array (Fig. 6b), while the higher frequency group demontsrated coherent resonances in the lower rows of the array (Fig. 6d). The group of resonances at the median frequencies could be characterized by a symmetry of DR

oscillations of two other SRRs.

Probe magnitude in A/m

resonance; and (e) the spectra of probe signals.

the right.

(e)

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 responses with respect to the central row i.e. the resonators located in upper and lower rows responded coherently, even though neighbouring columns demonstrated opposite phases of oscillations. So, it could be said that every other column responded coherently (Fig. 6c).

Fig. 5. Geometry, S21 spectra, and patterns of H-field magnitude (3d column) and phase (4th column) for SRR pairs arranged: (a) "face-to-face" and (b) "back-to-back ".

Fig. 6. (a) The spectra of the probe signals; inset shows the SRR array and the H-field probe locations; (b)-(d) typical patterns of the resonance oscillations in the array sampled at 6.64, 7.4 and 8.15 GHz (colour intensity shows the magnitude of resonance oscillations).

Despite of the above mentioned concerns about possible interaction between multi-element arrays and waveguide walls at the measurements we have performed comparison of simulation and measurement results for the transmission spectra at the placement of the 3 x 5 arrays of SRRs in the waveguide WR137. A good matching of the results presented in Fig. 7 is obvious, although both S21 spectra demostrate a smaller quantity of the transmission

Formation of Coherent Multi-Element Resonance States in Metamaterials 101

5.75 6.0 6.5 7.0 7.5 8.0 8.5 9.0 frequency, GHz

Fig. 8. (a) Geometry of a three- layer SRR array with marked probe locations; (b) coherent response of SRRs at 8.35 GHz; and (c) spectra of probe signals revealing enhanced splitting