Revisiting the wire medium: a resonant metalens
WWhy the wire medium with finite dimensions makes a resonant metalens
Fabrice Lemoult, ∗ Geoffroy Lerosey, and Mathias Fink
Institut Langevin, ESPCI ParisTech & CNRS, 10 rue Vauquelin, 75231 Paris Cedex 05, France (Dated: October 25, 2018)This article is the first one in a series of two dealing with the concept of ”resonant metalens” we recentlyintroduced [Phys. Rev. Lett. 104, 203901 (2010)]. Here, we focus on the physics of a medium with finitedimensions consisting on a square lattice of parallel conducting wires arranged on a sub-wavelength scale. Thismedium supports electromagnetic fields that vary much faster than the operating wavelength. We show that suchmodes are dispersive due to the finiteness of the medium. Their dispersion relation is established in a simpleway, a link with designer plasmons is made, and the canalization phenomenon is reinterpreted at the light of ourmodel. We explain how to take advantage of this dispersion in order to code sub-wavelength wave fields in time.Finally, we show that the resonant nature of the medium ensures an efficient coupling of these modes with freespace propagating waves and, thanks to the Purcell effect, with a source placed in the near field of the medium.
PACS numbers: 41.20.-q, 81.05.Xj, 78.67.Pt
I. INTRODUCTION
In a recent paper , we have introduced the concept of ”reso-nant metalens” which can be used to break the diffraction bar-rier. We defined it as a cluster of coupled resonators, arrangedon a subwavelength scale. When illuminated with broadbandelectromagnetic waves this macro-resonator radiates wavesthat contain information of an object placed in its near field.By measuring this radiation one can recover sub-wavelengthdetails of the object, and thus focus (or image) waves on adeep subwavelength scale, not limited by the Rayleigh crite-rion.In the original paper , we used metallic wires, half a wave-length long, as sub-wavelength resonators and we propose toexplain the physics of this example. We focus on the reso-nant metalens introduced initially consisting of a square lat-tice of N × N ( N = 20 ) identical metallic wires, alignedalong the vertical axis (defining the longitudinal z -direction).The length of the wires along z is equal to cm and theirdiameter to mm. The period lattice in both transverse direc-tions ( xy plane) is . cm, which is roughly equals to λ/ for the first intrinsic resonance of a single wire.The lens of interest here falls within the class of ”wiremedium”. We define the latter as a uniaxial medium formedby a periodic lattice of conducting wires with small radii com-pared to the lattice periods and the wavelength. In microwaveapplications this medium is often seen as an interesting arti-ficial dielectric, which presents the great advantage to have anegative effective permittivty . Some attention to wire me-dia has been paid in the realization of left-handed media ascomposite media made from lattices of long conducting wiresand split ring resonators . Usually, quasistatic models for thismedium are carried out for wave propagation perpendicularto the wires which show that, for electric-field polarizationalong the wires, the medium is characterized by a frequency-dependent effective dielectric constant.Belov and his colleagues has studied this kind of mediumin details, and especially the propgation along the wires. Theyshowed that this medium can be described by means of ananisotropic dispersive dielectric tensor. Also, they demon-strated that the wire medium supports propagating modes, so- called transmission line modes, which travel along the wires.During the whole article, we will focus only on these modeswhich explain our initial results and permit to give a simpleexplanation for our meaning.These modes which propagate in an adispersionless way inan infinite wire medium gives rise to dispersion when addingfiniteness along the wires. While in the original paper weinsisted on the fact that the dispersion is a key issue for sub-wavelength focusing/imaging with our method, we will give asimplified model for the physical mechanism responsible forit and a discussion justifying our simplification is presentedin appendix. We show why such a medium can be used tocode a sub-wavelength wavefield into a temporal/frequencialsignature thanks to its dispersion relation.This will lead us to note the link between our work andPendry’s designers plasmons , hence introducing the wavespropagating in a finite length wire medium as new spoof plas-mons. We will as well comment the very interesting effectwhich occurs at a specific frequency, known as canalizationregime , at the light of our simplified approach.The other key issue for a resonant metalens is the far fieldradiation of the sub-wavelength features. Usually, a sub-wavelength wave field is considered as evanescent, whichforbids its propagation toward the far field. But, the notionof evanescent waves, which is responsible for the diffractionlimit, is a mathematical formalism which only fits that of in-finite interfaces. Using finite transverse dimensions and end-ing on the medium presented initially , we show that the sub-diffraction details of the latter actually contributes to the far-field radiation.Eventually, it is commonly admitted that their contributionto the total radiation is much weaker than those of propagat-ing diffraction limited waves. Due to the resonant nature ofthe medium of interest, another physical mechanism whichis well known from opticians intervene: the Purcell effect .This effect increases the coupling of the source with radiatingmodes, and surprisingly it counterbalances the effect of theweak radiation of the sub-wavelength features. a r X i v : . [ phy s i c s . op ti c s ] N ov II. THE INFINITE WIRE MEDIUM: ADISPERSIONLESSPROPAGATION OF TEM BLOCH WAVES
We consider an infinite ”wire medium” consisting of infi-nite Perfect Electric Conductors (PEC) wires aligned alongthe z -axis, forming in the xy -plane a two dimensional squarelattice (Fig. 1). We assume a thin wire approximation, ie.their radii are really small compared to the lattice constant.We are looking for electromagnetic fields ( E ; H ) solutions ofthe following Maxwell problem and satisfying the boundaryconditions: (cid:126) ∇ × E = − µ ∂ H ∂t(cid:126) ∇ × H = ε ∂ E ∂t (1)where µ , ε are respectively the vacuum permeability andpermittivity. Here we have assumed that the electromagneticfields will propagate in vacuum only: the presence of metal-lic surfaces will be taken into account through boundary con-ditions only. Using the invariance of the medium along the z -axis and choosing a time dependence in e − iωt , we definetime-harmonic two-dimensional electric and magnetic fields E and H as: E ( x ⊥ , z, t ) = (cid:60) e (cid:16) E ( x ⊥ ) e i ( kz − ωt ) (cid:17) H ( x ⊥ , z, t ) = (cid:60) e (cid:16) H ( x ⊥ ) e i ( kz − ωt ) (cid:17) (2)where x ⊥ stands for the coordinates in the ( xy ) plane and k is the propagation constant of waves along the z -axis. E and H are complex valued fields depending on two variables(coordinates x and y ) but still having three components alongthe three axes. The medium of interest is an uniaxial mediumand we propose to decompose the electromagnetic fields withrespect to the geometry of the medium: E ( x ⊥ ) = E ⊥ ( x ⊥ ) + E (cid:107) ( x ⊥ ) e z H ( x ⊥ ) = H ⊥ ( x ⊥ ) + H (cid:107) ( x ⊥ ) e z (3)Eventually, due to the periodic nature of the medium inthe transverse plane, the problem becomes the well knownproblem in solid state physics which consists in looking forBloch waves solutions U k ⊥ (where U stands for E ⊥ , H ⊥ , E (cid:107) or H (cid:107) ) that have the form: U k ⊥ ( x ⊥ ) = e i k ⊥ . x ⊥ ˜ U ( x ⊥ ) (4)where k ⊥ is the so-called Bloch wave number (or quasi-momentum in solid state physics) and ˜ U ( x ⊥ ) is a functionthat satisfies the periodicity of the wire medium ( ˜ U ( x ⊥ ) =˜ U ( x ⊥ + R ) where R is a lattice vector). With these geomet-ric considerations the Maxwell system reduces to: k k c TETEM xz y k xy a~ /70 FIG. 1. (top) The geometry of the infinite wire medium: an in-finite square array of parallel perfectly conducting thin wires. Wefocus on the frequency range where the lattice parameter is sub-wavelength. (bottom) Dispersion curves for the TE and TEM typesof Bloch waves supported by the medium. The TEM Bloch modesare always propagating with a flat dispersion curve. For the geometryunder consideration ( r (cid:28) a (cid:28) λ ) the TM Bloch modes are alwaysevanescent . iωµ ˜ H (cid:107) e z = (cid:16) i k ⊥ + (cid:126) ∇ (cid:17) × ˜ E ⊥ iωµ ˜ H ⊥ = e z × (cid:16) ik ˜ E ⊥ − (cid:16) i k ⊥ + (cid:126) ∇ (cid:17) ˜ E (cid:107) (cid:17) − iωε ˜ E (cid:107) e z = (cid:16) i k ⊥ + (cid:126) ∇ (cid:17) × ˜ H ⊥ − iωε ˜ E ⊥ = e z × (cid:16) ik ˜ H ⊥ − (cid:16) i k ⊥ + (cid:126) ∇ (cid:17) ˜ H (cid:107) (cid:17) (5)Then, we study independently the three types of Blochmodes that can exist with respect to the wires: the transversemagnetic (TM) where ˜ H (cid:107) = 0 , the transverse electric (TE)where ˜ E (cid:107) = 0 , and the transverse electromagnetic (TEM)modes where ˜ H (cid:107) = ˜ E (cid:107) = 0 . The final equation to be solvedfor each type of mode becomes:TM : (cid:16) k ⊥ + k − (cid:0) ωc (cid:1) (cid:17) ˜ E (cid:107) − i ( k ⊥ .(cid:126) ∇ ) ˜ E (cid:107) − ∆ ˜ E (cid:107) = 0 TE : (cid:16) k ⊥ + k − (cid:0) ωc (cid:1) (cid:17) ˜ H (cid:107) − i ( k ⊥ .(cid:126) ∇ ) ˜ H (cid:107) − ∆ ˜ H (cid:107) = 0 TEM : (cid:0) k − ( ωc ) (cid:1) ˜ E ⊥ = (cid:126) (6)Solving the boundary value problem for TE and TEMmodes is easy since they do not have any electric field compo-nent along the wires, and a dispersion law is easily obtainablefrom the previous equations. For the TM modes this boundaryproblem is a tedious work which has already been performed .As a summary, the solutions of the Maxwell eigenproblem inthe first Brillouin zone must satisfy the following dispersionrelations (Fig. 1):TM : k ⊥ + k = ω c − k p TE : k ⊥ + k = ω c TEM : k = ω c (7)where k p is the so-called plasma wave number, a parameterthat depends on the lattice parameter a and the wires radii r and its expression can be found in the litterature .From now, we focus on the range of angular frequencies ω where the associated vacuum wavelength is really higher( > ) than the distance between two wires. In that caseand in the thin wire approximation, the plasma wavenumber isgreatly higher than the free space one and the TM waves arealways evanescent. The propagating TE modes along the z -direction are included inside the vacuum light cone. And theTEM Bloch modes are always propagating with a flat disper-sion curve independent of the norm of the Bloch wave number k ⊥ . Clearly, for the large k ⊥ corresponding to sub-wavelengthfeatures, the TE and TM modes are evanescent, while theTEM ones are propagating. We will only consider the propa-gating modes for this range, ie. the TEM Bloch modes. Thus,the fields inside the structure for a given k ⊥ can be describedby a unique scalar potential ˜ φ ( x ⊥ ) : E k ⊥ ( x ⊥ ) = − (cid:126) ∇ e i k ⊥ . x ⊥ ˜ φ ( x ⊥ ) H k ⊥ ( x ⊥ ) = 1 µ c e z × (cid:126) ∇ e i k ⊥ . x ⊥ ˜ φ ( x ⊥ ) (8) ˜ φ ( x ⊥ ) must satisfy the periodicity of the medium as well asthe boundary conditions at metallic surfaces. This potentialcan be seen as a spatial envelope that permits to verify theconstraints, and it multiplies a plane wave, as a consequenceof the Bloch theorem.During the rest of this article, we will always neglect thehigh diffraction orders due to this potential and consider it asa constant potential. We will always talk about plane wavesinstead of Bloch waves neglecting the Bloch potential. Oneimportant consequence of this hypothesis is that we can arbi-trarily change the coordinate axes without lack of generality,and especially we will often choose k ⊥ = k ⊥ e x .Another interesting result concerning the TEM Blochmodes is that the group velocity (cid:126)v g = (cid:126) ∇ k ( ω ) = c e z isconstant and oriented along the z -axis. Therefore any spatialwavepacket in a transverse plane propagates without distor-tion along the z -axis and one can use this property for guidingand imaging subwavelength features . III. THE FINITE WIRE MEDIUM ALONG THE WIRESAXIS: A DISPERSIVE MEDIUMA. Snell’s law for a wire medium-air interface
Unfortunately, the infinitely extended wire medium cannotexist and introducing a boundary condition leads to interestingresults. Here, we consider a TEM Bloch wave with transversewave number k ⊥ = k ⊥ e x traveling in the half space z < L/ filled out with a wire medium. This Bloch wave encountersa boundary with vacuum at z = L/ . The incident angle θ i is defined by k ⊥ /k . This wave generates a reflected wavewith angle θ r and a transmitted wave with angle θ t . The firststep consists in finding an equivalent of the Snell’s law forthis interface. Considering the conservation of the tangentialcomponent of the wave number at the interface (here k ⊥ ) andintroducing the dispersion relation on both sides of the equal-ity we obtain: θ r = − θ i tan θ i = sin θ t (9)From this refraction’s law the critical incident angle θ c = π/ appears: when the incident angle is higher than θ c (ie. forevanescent waves in the air region) a total internal reflectionis observable.In order to simplify the calculation for continuity of tan-gential components, as explained in the previous part, we willneglect the evanescent TM wave that is generated in the wiremedium at the reflection and the high diffraction orders. Adiscussion concerning those approximations is presented inappendix. Then, according to equation (8), the electric fieldson both sides of the interface write: E ( i ) = E e i ( k ⊥ x + k ( z − L/ E ( r ) = rE e i ( k ⊥ x − k ( z − L/ (10) E ( t ) = tE e i ( k ⊥ x + κ ( z − L/ cos θ t − sin θ t And the associated magnetic fields write: H ( i ) = E µ c e i ( k ⊥ x + k ( z − L/ H ( r ) = rE µ c e i ( k ⊥ x − k ( z − L/ − H ( t ) = tE µ c e i ( k ⊥ x + κ ( z − L/ (11)where k = ωc stands for the free space wave number, and κ = √ k − k ⊥ if θ i ≤ θ c or κ = i √ k ⊥ − k if the total internalreflection occurs. Applying the continuity of the tangentialcomponents of the fields at the interface and introducing therefraction law of equation (9) yields to the following reflectioncoefficient: r = − k − κk + κ (12)Interestingly, in the peculiar case of total internal reflection, r becomes a complex coefficient with unity norm, and the re-flection induces a phase shift: ϕ ( k ⊥ , k ) = π − (cid:18) √ k ⊥ − k k (cid:19) (13)This phase shift is a well-known phenomenon for total internalreflection and one of its consequences is the Goos-H¨anchenshift along the x -direction. Here, this phase shift has animportance when adding a second boundary at z = − L/ .Now, the wire medium with finite dimension consists of wiresaligned along the z -direction with the same length L (Fig. 2).With this new constraint the field inside the medium experi-ences a Fabry-Perot like resonance, or equivalently a guidedpropagation along x -axis, when the optical path length δ for around trip satisfies: δ = 2 kL + 2 ϕ ( k ⊥ , k ) = 2 nπ (14)where n is an integer. With equation (13) and (14) we canwrite a dispersion relation for the guided waves inside thebounded wire medium for the first resonance ( n =1): tan (cid:18) k L (cid:19) = (cid:112) k ⊥ − k k (15)This implicit theoretical dispersion relation is plotted interms of the first Fabry-Perot frequency f m versus k ⊥ (bothnormalized to f and k corresponding to an ideal Fabry-Perotof length λ/ ) in Figure 2. The asymptotes of the light line atlow frequencies and the intrinsic resonance frequency at largevalues of k ⊥ are shown.Rigorously, the TM evanescent waves as well as the highdiffraction orders (and thus the TE diffracted waves inside thewire medium) induced at the reflection adds small correctionsto the phase shift. This would lead to a more complicated dis-persion law that depends on the wire medium parameters. Forsimplicity and because this dispersion is in a sufficiently goodagreement with simulation results (Fig. 7), we have preferrednot to introduce those waves. More sophisticated dispersionlaws are presented in the appendix. We believe our calcula-tion that only takes into account the evanescent field at reflec-tion is the easiest way to understand why the z -bounded wiremedium introduces the dispersion. Because this dispersion isa key issue for the concept of resonant metalens, we have cho-sen to make those approximations for the scope of this article,but obviously one can add small corrections to this law. B. Fields profiles inside the medium
In order to describe precisely these propagating waves in-side the wire medium we consider that the region between z = ± L/ is filled out with a wire medium (Fig. 2). Weknow that the field inside this region is TEM with respect tothe wire but the invariance along z is broken. Still consideringthe first order of diffraction, neglecting TM evanescent waves,and choosing the x -axis in order to have k ⊥ = k ⊥ e x , we canexpress the fields inside the wire medium as: E TEM = e i k ⊥ . x ⊥ E ( z )00 H TEM = 1 iµ ω e i k ⊥ . x ⊥ d E / d z (16)And, from the propagation equation (5) of the TEM Blochmodes inside the wire medium it follows: E ( z ) = C sin (cid:16) ωzc (cid:17) + C cos (cid:16) ωzc (cid:17) (17)In the air region on both sides of the wire medium the elec-tromagnetic fields decompose onto a plane wave that we canwrite using the previous notations: E vac = E , e i ( k ⊥ . x ⊥ ± κ ( z ∓ L )) ± κ/k − k ⊥ /k H vac = E , iµ c e i ( k ⊥ . x ⊥ ± κ ( z ∓ L )) (18)Then applying the continuity of the tangential fields at theinterface we reach the matrix form equation: sin ( kL/
2) cos ( kL/ − κ/k − cos ( kL/
2) sin ( kL/ − i − sin ( kL/
2) cos ( kL/
2) 0 κ/k − cos ( kL/ − sin ( kL/
2) 0 − i C C E E = (cid:126) (19)The system has a non-zero solution only if its matrix de-terminant is null. The first solution of this cancellation givesanother expression of the dispersion relation of equation (15): cos (cid:18) k L (cid:19) = kk ⊥ (20)This expression is actually the same as the one presentedpreviously in equation (15), and one can develop the sine func-tion from its cosine and recover the tangent function. Fromthis expression, we can eventually express all of the electro-magnetic components. In the wire medium region it gives: k / k f m / f L eff L /2 Goos-Hänchen shift E H E z xz y FIG. 2. (top) The geometry of the infinite wire medium with finitelength along the wires. (bottom) The dispersion relation of the firstresonance frequencies f m (normalized to the original Fabry-Perotfrequency f ) versus the transverse Bloch wave number k ⊥ (nor-malized to the original Fabry-Perot wave number k ) from equation(15),(20) or (25). The asymptotes of the light cone and the originalresonance frequency are shown (dashed-lines). The inset shows howthe three ways of calculation (reflection coefficient, guided wavesexpression and effective Fabry-Perot length) leads to the same result. E TEM = E sin ( kz ) e i ( k ⊥ . x ⊥ ) e x H TEM = − iE µ c cos ( kz ) e i ( k ⊥ . x ⊥ ) e y (21)One can notice that for the first resonance the electric fieldis maximum at both ends of the wire medium, while the mag-netic field is maximum at half the wires length, and is nearzero at both ends (fig. 2). And in the vacuum regions for highBloch wave numbers (ie. the waves are evanescent in air) itgives: E vac = E e i k ⊥ . x ⊥ e ∓ √ k ⊥ − k ( z ∓ L ) ± (cid:112) − k /k ⊥ i H vac = iE µ c kk ⊥ e i k ⊥ . x ⊥ e ∓ √ k ⊥ − k ( z ∓ L ) e y (22) C. Energy density considerations
However, we can use a more intuitive approach to deter-mine the dispersion law. In fact, depending on the trans-verse wavevector k ⊥ , the Bloch waves penetrate in air up toa distance whose inverse is the quantity previously introduced (cid:112) k ⊥ − k . Again, a more sophisticated model would addcorrection to this penetration depth taking into account highdiffraction orders and TM waves (see appendix). Each propa-gating mode in the structure experiences the Fabry-Perot res-onance, and the length of the cavity is not simply the length L of the wire medium, because of the penetration of the modesinside air. A common way to describe the effective length ofthe real cavity is to calculate the quantity: L eff = (cid:20)(cid:90) ∞−∞ u ( x, y, z ) d z (cid:21) (cid:90) ∞−∞ [ u ( x, y, z )] d z (23)where u ( x, y, z ) is the energy density of a given mode. Us-ing the fact that the modes are TEM inside the structure, thedensity of energy inside the structure is constant along the z -direction, from − L/ to L/ . For the sake of simplic-ity, we set this value to , which does not hamper our de-scription since this value simplifies in L eff . The energy den-sity decreases exponentially at the interfaces z = − L/ and z = L/ and taking into account the three non-zero compo-nents of the fields, the calculus simply gives: L eff = L + 2 (cid:112) k ⊥ − k (24)This expression has merit to show explicitly the increaseof the Fabry-Perot cavity length due to penetration depths atinterfaces. We finally consider that the resonance condition isobtained if L eff = nλ/ for the n th Fabry-Perot like condition,which leads after rearrangement to the dispersion relation ini-tially presented for the first Fabry-Perot resonance: f f = 1 + 2 π (cid:34)(cid:18) k ⊥ k (cid:19) − (cid:18) ff (cid:19) (cid:35) − (25)where k and f stands for the original first Fabry-Perot likeresonance, namely k = π/L and f = c/ L . D. Link with designer plasmons
The dispersion curve deserves some comments. Indeed,this relation is similar to the ones characterizing surface plas-mons polaritons (SPPs) and we have demonstrated the exis-tence of guided modes inside the wire medium (with finitelength) with high wavenumbers. But, at low frequencies,where the limit of perfect conductor is applicable, it is wellknown that perfect metals do not support electromagnetic sur-face modes, forbidding the existence of SPPs.However, it has been demonstrated theoretically andexperimentally that bound electromagnetic surface wavesmimicking SPPs can be sustained even by a perfect conductor,provided that its surface is periodically corrugated. Two dis-tinct geometries of modulation of the flat perfectly conduct-ing surface have been introduced : a one-dimensional arrayof grooves and a two-dimensional hole array. The frequen-cies of the supported modes scale with the geometrical sizeof the corrugations in the perfect conductor approximation.Maier and his colleagues also demonstrated the existence ofpseudo surface waves in the low frequency regime: structuresconsisting of closely spaced metal rods in one dimension canguide electromagnetic waves. In the present case, the sameconfiguration has been investigated but in the two dimensionalcase, and we have demonstrated the same guiding properties.Interestingly, the equation (15) is actually the same disper-sion law as the one obtained for the spoof plasmons . Again,in both cases the high diffraction orders and the mode conver-sions are neglected when calculating the dispersion law, andthis is the reason why the two laws are the same. Because ofthe matching of the dispersion laws at the first order, one canimagine to couple the two kinds of structure. Similarly, thethickness of the layer that supports these plasmon-like wavesscales with the operating frequencies since the height of thewires is near λ/ . But, here, only the length of the wires hasan influence on the resonance frequency. It is not the case forthe array of grooves where the groove’s width and the distancebetween two of them have an impact on the dispersion curve.In the array of pinholes, the fundamental mode inside a holeis considered which implies a dependence of the hole dimen-sions in the dispersion. Here, the distance between wires aswell as wires diameters do not influence the dispersion lawand the guiding properties of the medium do not seem to besignificantly altered when wires are displaced from their orig-inal positions. Another interesting property is that the wiremedium is essentially made out of vacuum, thus losses in-duced by metallic surfaces are less important than in the cor-rugated surfaces or in the surface with pinholes. E. Link with canalization regime
Eventually, we would like to add comments concerning theresults obtained by Belov and his colleagues . They haveused the wire medium for performing a pixel to pixel imag-ing with a resolution thinner than the diffraction limit. Thesame near field lens has also recently reached the near infra-red spectrum and sub-wavelength imaging has been provedexperimentally using gold nano-wires. In the two cases, theapproach is monochromatic and the operating frequency oftheir lens corresponds to the intrinsic resonance of a singlewire, namely f . At this exact frequency the field distributionat an interface is recovered at the other interface, and by scan-ning near field of one interface they perform sub-wavelengthimaging. We would like to prove that their results correspond-ing to the canalization regime are consistent with the ones de- veloped here. From previous calculations, introducing equa-tion (15) into equation (12) leads to the following reflectioncoefficient: r ( f ) = − − i tan( πf L/c )1 + i tan( πf L/c ) (26)Therefore, at the intrinsic resonant frequency of a singlewire f this reflection coefficient is equal to − and does notdepend on the Bloch wavenumber k ⊥ since the dispersion re-lation has no solution at this frequency. Because this coeffi-cient is real, the reflection does not imply any phase shift andespecially the Goos-H¨anchen shift is null. As we have seenbefore, the group velocity for the TEM Bloch waves inside thewire medium is oriented along the wire axis, so if there is noGoos-H¨anchen shift at the reflection, there is no propagationin the transverse plane. At f , whatever the Bloch wavenum-ber k ⊥ considered the reflected wave propagates along thewires as well as the incident one, hence giving rise to the socalled canalization regime. This regime implies that the elec-tric field has the same profile at both interfaces. When placinga sub-wavelength source at an interface they managed to im-age it half a wavelength away at the other interface with aresolution far better than the diffraction limit . F. Coding of sub-wavelength details
The dispersion is a key issue for the concept of resonantmetalens we introduced . The importance of this dispersionmanifests when a sub-wavelength source is placed near oneof the two interfaces. Actually, the two interfaces correspondto the object plane of the lens and the lens permits to obtainthe field profile inside one of these transverse planes. Thisprofile decomposes with a unique set of amplitude and phaseon the modes supported by the wire medium. Typically, it canbe expressed in terms on k ⊥ since it corresponds to a Fouriertransform: E ( x ⊥ ) = (cid:90) (cid:90) ˜ E ( k ⊥ ) e i k ⊥ . x ⊥ d k ⊥ (27)where in the present case the integration domain is limited tothe first Brillouin zone of the lattice, that is to say when oneof the component of the Bloch wavenumber reaches π/a .This formalism makes sense when considering broadbandilluminations. Interestingly, from previous calculations andignoring the degeneracy, a given Bloch wavenumber k ⊥ isassociated to its own angular frequency ω . For the vacuumdispersion relations this is also the case but here k ⊥ can be as-sociated to a spatial scale that is deeply sub-wavelength com-pared to the associated vacuum wavelength. This propertypermits to transform the function ˜ E ( k ⊥ ) into a function thatdepends only on frequency: this is the frequency coding ofthe sub-wavelength details. This property means that, whenthe medium is illuminated from one of the two object planes,registering the frequency spectrum of a wave packet gives thespatial Fourier transform of the source. The spatial features ofthe source can now be expressed in terms of ω : ˜ E ( k ⊥ ) = F ( ω ) (28)Eventually, it is well known that time and frequency are twodistinct representations of the same phenomena. A broadbandillumination means the emission of a short pulse, and the fre-quency signature F ( ω ) for the sub-wavelength features meansa temporal one f ( t ) . The imaging procedure associated to thisconcept is the scope of the next article, but at this time one caneasily understand that the resolution must be limited by thedistance between two wires. By exploiting the dispersion ofresonant modes with broadband illuminations we can recordthe sub-wavelength details of an object. In the next part wewill see that this information can be carried toward far field. IV. A WIRE MEDIUM WITH FINITE TRANSVERSEDIMENSIONS: FAR FIELD RADIATION OFSUB-WAVELENGTH FEATURES
Until now, everything happened inside or in the near fieldof the wire medium. Now, we consider the previous wiremedium but with finite dimensions in both transverse direc-tions (Fig. 3). Introducing this finiteness has a first conse-quence which is the quantization of the transverse wavenum-bers. In solid state physics, the quantization is done by intro-ducing BornVon Karman boundary conditions which con-sists in cycling a crystal. When studying crystals the numberof atoms considered is large enough to allow this approxima-tion. In the present case, the number of wires in both trans-verse directions is only and this boundary condition is notwell suited.But, introducing finiteness results in the same quantizationeffect. At each interface, the plasmon-like wave traveling in-side the wire medium encounters a boundary with vacuumand gives birth to a reflected one. The analytical expressionof the reflection coefficient is not easy to calculate since thefield in the air region is the superposition of many plane wavesdue to the finiteness along the wires. Then, the interferencesinside the wire medium are constructive only for a discretespectrum of Bloch wavenumbers which characterize the sta-tionary eigenmodes of the system. Thus, due to the finite-ness of the array of wires, the transverse Bloch wave num-bers are quantified: k ⊥ = πD ( m. e x + n. e y ) , with integers ( m, n ) ∈ [[0; N − and D the size of the medium in the x and y directions, D = a ( N − . Again, from the previ-ous section, it is important to keep in mind that each k ⊥ ( m,n ) is associated to its own resonant frequency f ( m,n ) , ensuring afrequency (or equivalently a time) coding of spatial features.The second consequence of these interfaces is a conver-sion from the plasmon-like modes to propagating waves inthe free space region. Intuitively, with an approach of coupledresonators, one easily understands that each wire generates amonopolar radiation (for the first resonance) which is essen-tially z -polarized. The wire medium can therefore give rise toa superposition of monopolar radiations. With the approach L /2 xz y N wires N wires EH FIG. 3. The geometry of the wire medium with finite dimensionsin the 3 spatial directions. Due to the finiteness this medium radiatesfield toward the far field and superimposed on the medium we haverepresented a sketch of a radiation pattern for the first resonance.This is just an example since the next figure shows that the directivitypattern in the ( xy ) plane can be different. we used until now, the simplest way to explain the radiationcomes from the fields profiles inside the wire medium (insetof figure 2). We have already seen that the fields cannot exitthe structure from the top or the bottom interfaces becausethe waves are evanescent. The only way to escape the wiremedium is from lateral interfaces. Since the x and y com-ponents of the electric field present an odd profile for the firstFabry-Perot resonance, they cannot contribute to far field radi-ation. While the z -component of the electric field can be mod-eled by two small dipoles (cid:126)p , positioned at z = ± L/ . Thefield radiated by the superposition of these two infinitesimaldipoles writes, in the dipolar approximation and in sphericalcoordinates: E = µ ω p sin Φ4 πr e ikr (cid:16) e i kL sin Φ + e i ( − kL sin Φ+ ϕ ) (cid:17) e Φ (29)where Φ corresponds to the zenith angle, and ϕ stands forthe phase difference between the two dipoles. Depending onthe parity of the Fabry-Perot resonance order the two dipolesare in phase or out of phase ( ϕ = 0 or π ). Typically, for thefirst Fabry-Perot resonance, or for the field profile presentedin figure 2, the two dipoles are in phase, ie. ϕ = 0 . Thus,for the first Fabry-Perot resonance, the radiated amplitude ismaximum in the plane z = 0 and it can be studied as a twodimensional problem.Now we consider that we have such small dipoles on bothinterfaces of the wire medium with a 2D spatial distributiongiven by the z -component of the TEM Bloch mode at the in-terfaces. We perform the projection of this mode onto thebasis of 2D free space plane waves. The E z field at the inter-faces of the structure has the periodicity given by the Blochwavenumber and has the form: E k ⊥ z ∝ sin (cid:18) k ⊥ x ( x − D (cid:19) sin (cid:18) k ⊥ y ( y − D (cid:19) Π D ( x, y ) (30)where Π D is the 2D rectangle function of width D . −20 0 20−20020 −20 0 20−20 0 20−20020−20020 −1 0 1−101 −1 0 1−101 −1 0 1−101 k . e y / k k . e y / k k .e x /k k .e x /k k . e y / k FIG. 4. The two dimensional projection of three distinct Blocheigenmodes onto the free space waves. The left colormaps show thatthe finiteness of the wire medium in the transverse plane lead to asuperposition of 4 sinc functions. The right ones is the projectionof the left ones onto the light cone: when considering object withfinite dimensions a certain amount of energy can reach the far field.These colormaps also show that this representation gives informationabout the directivity pattern of an eigenmode. Depending on its k ⊥ monopolar, dipolar or quadrupolar radiation can occur. In the Fourier domain, this expression transforms into asum of four sinc functions positioned at the corners of a rect-angle of width k ⊥ x and length k ⊥ y (Fig. 4). Such a repre-sentation contains all information about the propagating fieldsgenerated by the sub-wavelength modes. The intersection ofthe latter with a disk of radius ω/c (green line in figure 4)gives the efficiency of the radiation as well as the directivitypattern.Averaging this intersection along the disk resides in an es-timation of the radiated energy in the far field for a givenmode. The full calculation of the radiated energy residesin a projection onto a set of Hankel functions and is quitehard to perform, but an estimation of this quantity can bedone. Due to the fact that we sum 4 sinc functions whoseoperands are proportional to Dk ⊥ x,y , and due to the fact that the transverse dimension is lower than the operatingwavelength, the radiation’s efficiency in energy is roughlyequal to / (16 π ) . ( DLk ⊥ x k ⊥ y ) : the smaller the transversewavenumber the higher the radiated amplitude. Fortunately,we will see later that the Purcell effect compensates for thislow radiation efficiency.In the same time, the spectral representation of the modesgives information about the directivity pattern of the radiatedfield. The intersection between this spectrum and the circleof radius k directly gives the radiation pattern of the modes.A point on the circle with coordinates ( cos θ , sin θ ) permitsto obtain the radiated amplitude of the mode in the directiondefined by the angle θ . The directivity pattern is thereforefully described by this representation. It is important to keepin mind that this representation is in the ( xy ) -plane and thatthe radiation in a plane containing ( Oz ) is like the one of asingle dipole antenna. To underline this aspect we representthe projection of the previous k -space spectra onto this circle(Fig. 4).Depending on the mode parity, or equivalently of the parityof m and n , the Fourier spectrum has to keep the same parity.For a given direction in the transverse plane, if the integerdescribing the wave number component is odd, the spectrumis even: it results in a null a value for the spectrum at theorigin. On the contrary, an even integer does not impose anull value at the origin, and it does not cancel over because inthe present case the transverse dimension D is lower than halfa wavelength. Taking into account the two directions it resultsin four directivity pattern possibilities: • the odd-odd modes give monopolar radiation; • the odd-even (resp. even-odd) modes give a dipolaralong x (resp. y ) radiation; • the even-even modes give a quadrupolar radiation.Therefore the resonant metalens possesses four spatial de-grees of freedom, representing as many information channelsexploitable for imaging and focusing . At a frequency wheresome modes are degenerate, registering far field in differentdirections allows a discrimination between modes, by the wayincreasing the number of subwavelength information regis-tered.This type of directivity patterns can also be interpretedin terms of corner radiation which have been introducedin acoustics. In calculating radiation from vibrating squareplates, a specific kind of modes, sometimes termed cornermodes for the reason that only the corner quarter-cells con-tribute significantly to the radiation, have been described.Those modes satisfy the condition k < mπ/D and k The idea we introduced in the original paper was to usethe radiation of these subwavelength features to image an ob-ject placed near the resonant medium. Thus the last issue tobe resolved is the coupling between the object and the lensitself. Spontaneous emission is not only an intrinsic propertyof the atoms or molecules in a specific material but is alsogoverned by the electromagnetic properties of the host. In the1950s, E.M. Purcell demonstrated that the rate of sponta-neous emission for a light source depends on its environment:it can be increased when placing the source in a resonant cav-ity. This phenomenon is nowadays well known as ”Purcelleffect”.Later, E. Yablonovitch predicted theoretically that a ”pho-tonic” material can control the rate of radiative recombinationof an embedded light source. This led to the huge domain ofphotonic bandgap materials where this effect is used for in-hibiting the spontaneous emission of light.Also, many theoretical works concern the effect of an ab-sorbing dielectric on spontaneous emission and level shifts ofan embedded atom, using a Green’s function approach , orusing a macroscopic Lorentz’s cavity model , or using equa-tions of motion for a collection of two-level atoms embeddedin a dielectric host medium . As an atom in a dielectric issurrounded by many other two level resonators, these worksare similar to the concept of resonant metalens we study here.Interestingly, we have already discussed the frequency shiftin the case of the resonant metalens and here we discuss theincrease of the spontaneous emission rate. Namely, in the res-onant metalens case, the Purcell effect implies an increase ofthe coupling between the source and the wire medium sinceit acts as a cavity surrounding the source and modifies the so-called local density of states (LDOS).The computation of the eigenmodes readily yields the qual-ity factor of the modes resonance, a factor influencing the 200 220 240 260 280 300 320 340 360 Frequency (MHz) R e t u r n l o ss 200 220 240 260 280 300 320 340 3600.60.70.80.911.1 Frequency (MHz) S p a r a m e t e r ( u . a . ) FIG. 6. Measured return loss for a small dipole with (blue) andwithout (red) the wire medium. Due to the resonant nature of thewire medium a completely reactive dipole becomes radiative, andthis efficiency grows with the resonance Q factor, and thus with fre-quency here. The discrete nature of the blue curve still demonstratesthat the subwavelength features are coded in the field spectrum. spontaneous emission. From the previous calculation, andin the deeply sub-wavelength limit, we showed that the ra-diation’s efficiency is roughly proportional to ( DLk ⊥ x k ⊥ y ) :the lower the transverse wave number, the higher its radiatedenergy. In other words the resonance quality factor, which isthe ratio between the energy stored in the resonant mediumand the radiated energy, grows as ( DLk ⊥ x k ⊥ y ) . As statedby E.M. Purcell , it results in a better rate of spontaneousemissions for the highest transverse wave numbers.This effect manifests itself by a better coupling betweenthe source in front of the lens and especially for the most sub-wavelength modes. The resonant nature of the modes matchesthe impedance of the small electric dipole, or equivalently, thePurcell effect compensates for the weak radiation of the deepsub-wavelength modes. As a summary, we have seen that thePurcell effect implies a rate of spontaneous emission propor-tional to the quantity ( Dk ⊥ ) , while the radiation efficiencyof a given mode is proportional to the inverse of this quantity.Overall, the two effects counterbalance themselves ensuringthat all of the subwavelength details radiate the same amountof energy. And the higher the Bloch wavenumbers, the thinnerthe bandwidth carrying the energy.On a temporal point of view, the Q factor is the oppositeof the attenuation time of the system: a higher quality factorimplies a lower damping, and so high Q modes oscillate forlonger times. This is why in the original paper of metalenses we presented a curve of the lifetimes of the modes inside thestructure versus their Bloch wavenumbers k ⊥ . Here, we haveexplained why the data were growing with k ⊥ . The linear fitinitially obtained results from the short time accessible for thetime domain simulation compared to the extremely high Q factor attainable.In order to seal the validity of this physical description, wemeasured the return loss of a small dipole with a network an-alyzer (Agilent N5230C) with and without the presence of thelens (Fig. 6). The return loss is the ratio of the power of theoutgoing signal to the power of the signal reflected back fromthe small dipole. If the return loss is equal to 1, the wholeinput power is reflected back meaning that there is no radi-ated field. While when this quantity decreases, some amountof energy has been transmitted by the dipole: this energy0corresponds to the radiated energy. Without the lens, it re-sults in a flat curve in amplitude demonstrating that the smalldipole is totally reactive, that is to say, it does not radiate anyfield. In the presence of the resonant metalens, the curve de-creases down to a value of 0.7, meaning that the source isbetter matched due to the Purcell effect as stated before. Thebetter coupling between the source and the resonant mediumpermits to compensate the less efficient radiation of the mostsubwavelength modes: the two effects counterbalance them-selves and each mode eventually radiate the same amount ofenergy. VI. CONCLUSION In this article, we have explained the physical mechanismof our original paper concerning the concept of resonantmetalens . The intuitive property of resonance splitting whenconsidering coupled resonators has been investigated througha wave approach. Starting from the wavefield inside an infi-nite wire medium, we have seen that a wire medium with fi-nite dimensions can be seen as a macro resonator that containsmany resonance frequencies. Interestingly, each of these reso-nances corresponds to an eigen wavefield that fluctuates on itsown sub-wavelength scale. Then, we have shown how thesemodes can radiate into the far field region information that isusually contained in the evanescent spectrum. Eventually, byplacing such a resonant medium near a sub-wavelength sourcewe have seen that the Purcell effect permits an efficient radia-tion.These physical mechanisms permit to get information ofa source on a smaller scale than the wavelength, stored as atime/frequency signature, in the far field region. In the nextarticle we will explain how one can use these properties forbreaking the diffraction barrier in a focusing/imaging schemeand study the main limitations as well as losses impact. Appendix: Sophisticated dispersion law When calculating the dispersion law of guided modes in-side the wire medium, we have took sides for neglecting boththe TM evanescent waves and the high orders of diffractiongenerated at the reflection. For the range of frequency con-sidered, the TM waves are strongly evanescent because theplasma wave number k p is really higher than any other wavenumber considered . And the high order of diffractions arealso strongly evanescent because the fundamental one is al-ready evanescent. All of these waves should modify the re-flection coefficient of equation (12) by an additional factorwhich affects phase, and thus dispersion.A dispersion law of the guided modes considering the TMwaves is presented in the appendix of ref. . To obtain this dis-persion law the authors used an effective medium approachwhich leads to an effective anisotropic permittivity along thewire axis . It permits to obtain the Bloch wave number ofthe TM modes: κ TM = (cid:112) k − k ⊥ − k p , which is a com-plex quantity for the range of frequency considered here. To keep the formalism of their article, it is convenient to usethe notation κ TM = iγ TM . Also, by analogy, we denote thequantity κ = √ k − k ⊥ = iγ , where γ is a real quan-tity for the waves we consider here (evanescent waves in freespace). Then, an Additional Boundary Condition is intro-duced to solve the boundary problem at the air/wire mediuminterface , and the law obtained can be summarized in: tan (cid:18) k L (cid:19) = γ app k (A.1)with, γ app = γ k ⊥ + k p k p + γ TM k ⊥ k p tanh( γ TM L (A.2)This dispersion relation is actually the same as the one pre-sented in the manuscript with a correction due to the TMwaves which appears in γ app . In the limit of thin wires andfor frequencies well below the cutoff of TM waves (which isthe case considered here), the plasma wave number k p is ex-tremely higher than any other wave number considered. Inthat limit case we find that γ app = γ and the dispersion lawpresented in equation (15) remains true. With the geometricalparameters we took in this example, the difference betweenthe two laws is less than 1 % (Fig. 7): this is the reason whywe chose not to take into account the strongly evanescent TMwaves.As seen in the expression of γ app which is always greaterthan γ , the TM waves tend to increase the apparent wavenumber associated to the evanescent decrease in air, or equiva-lently to decrease the penetration depth inside air. The guidedmodes are much more weakly confined to the interface thanthe case without considering TM waves. As a consequence,for a given transverse wavenumber k ⊥ it results in a higherfrequency for the guided mode than the dispersion law pre-sented in the paper.We would like to argue that the dispersion law obtainedwith 3D time domain simulations (Fig. 7) gives a distinctcomportment for large k ⊥ than the one expected when addingthe TM waves. Indeed, while the correction due to the TMwaves tend to increase the frequency, the simulation resultsshow that a decrease is expected. Actually, the periodic na-ture of the wavefield and thus the Bragg waves have beenneglected, a remaining problem with effective medium ap-proach. Here we propose to do the calculation of the disper-sion law taking into account some high orders of diffractionwith a method inspired from ref. .We consider a TEM Bloch wave with transverse wave num-ber k ⊥ = k ⊥ e x traveling in the half space z < filled outwith a wire medium. This Bloch wave encounters a boundarywith vacuum at z = 0 and gives rise to a reflected TEM Blochwave inside the wire medium region. The superposition of the2 TEM Bloch waves writes at z = 0 :1 ⊥ /k f m / f only TEMTEM + TMTEM + diffraction3D simulation FIG. 7. The dispersion laws of the first resonance frequencies f m (normalized to the original Fabry-Perot frequency f ) versusthe transverse Bloch wave number k ⊥ (normalized to the originalFabry-Perot wave number k ) from equation (15) only consideringTEM waves in the wire medium (blue solid line), from ref. withour parameters considering TM and TEM waves (red dashed line),and from equations (A.1) and (A.9) considering only TEM wavesbut high diffraction orders (green dashed dotted line). Superimposedwith stars is the results obtained from 3D simulations. E TEM = (1 + r ) e ik ⊥ x F ( x, y )00 H TEM = 1 µ c (1 − r ) e ik ⊥ x F ( x, y )0 (A.3)The function F ( x, y ) stands for the Bloch envelope of theTEM wave inside the wire medium, and notably this functionis null over the wire section of a unit cell. One can notice thatthe y (respectively x ) component of the electric (respectivelymagnetic) component of the TEM wavefield has been ne-glected: the diffraction orders along the y -axis are neglected.This implies that the reflected and transmitted TE modes willbe neglected too, but in the general case they must be takeninto account. We will also neglect the TM reflected wavesbecause their influence has been discussed earlier. Thus, thetransmitted field in the vacuum region, with these approxima-tions, writes at z = 0 : E vac = ∞ (cid:88) n = −∞ t ( n ) e ik ( n ) ⊥ x k ( n ) z /k − k ( n ) ⊥ /k H vac = 1 c ∞ (cid:88) n = −∞ t ( n ) µ c e ik ( n ) ⊥ x (A.4)with k ( n ) ⊥ = k ⊥ + nπa and in the range of frequency consid-ered, we have again k ( n ) z = iγ n = i (cid:113) k ( n ) ⊥ − k . To eliminate the unknowns, we use the fact that the tan-gential components of the electric field must be continuousat the vacuum/wire medium interface over the entire unit cell(which represents the area Ω ), while the magnetic field com-ponents are continuous only at the air regions of the unit cell(the area Ω minus the area D corresponding to the wire sec-tion). Matching the x component of the electric field, multi-plying by e − ik ( m ) ⊥ x , and integrating over the area Ω yields: t ( m ) = kk z ( m ) (1 + r ) I m (A.5)where, I m = 1 a (cid:90) (cid:90) Ω F ( x, y ) e ik ⊥ x e − ik ( m ) ⊥ x d S (A.6)Again, matching the y component of the magnetic field,multiplying by F ( x, y ) e − ik ⊥ x , and integrating over the area Ω − D gives: ∞ (cid:88) n = −∞ t ( n ) I n ∗ = (1 − r ) J (A.7)where, J = 1 a (cid:90) (cid:90) Ω −D F ( x, y ) d S (A.8)Then we can extract the reflection coefficient of the TEMwave, and as done in the manuscript when adding a secondalong the wires, it eventually gives the same dispersion law asequation (A.1) with in this case: γ app = J ∞ (cid:80) n = −∞ | I n | γ n (A.9)To evaluate this quantity, we need to give an expression forthe function F ( x, y ) . For simplicity we have chosen that thisfunction is equal to in the domain ω − D , and null in the do-main D . The dispersion law obtained considering the 2 firstdiffraction orders ( +1 and − ) is presented in figure 7. Onecan notice that this dispersion law gives a better fit to the re-sults obtained from the full 3D simulations. Obviously, thiscalculation is not an exact one but it has the merit to showthat taking into account the diffraction orders tends to increasethe apparent penetration depth inside air, thus giving a bettercomportment for the high transverse wave numbers. Never-theless, this calculation greatly complexifies the understand-ing for giving a correction less than 1 %: this is reason whywe also neglect the high diffraction orders during the wholemanuscript.2 ∗ [email protected] F. Lemoult, G. Lerosey, J. de Rosny, and M. Fink, Phys. Rev. Lett. , 203901 (2010) J. B. Pendry, A. J. Holden, W. J. Stewart, and I. Youngs, Phys.Rev. Lett. , 4773 (Jun 1996) R. A. Shelby, D. R. Smith, and S. Schultz, Science , 77 (2001) P. Belov, S. Tretyakov, and A. Viitanen, Journal of Electromag-netic Waves and Applications , 1153 (2002) P. A. Belov, R. Marqu´es, S. I. Maslovski, I. S. Nefedov, M. Sil-veirinha, C. R. Simovski, and S. A. Tretyakov, Phys. Rev. B ,113103 (Mar 2003) P. A. Belov and M. G. Silveirinha, Phys. Rev. E , 056607 (May2006) J. B. Pendry, L. Martin-Moreno, and F. J. Garcia-Vidal, Science , 847 (2004) P. A. Belov, Y. Hao, and S. Sudhakaran, Phys. Rev. B , 33108(2006) P. A. Belov, C. R. Simovski, and P. Ikonen, Phys. Rev. B ,193105 (May 2005) E. Purcell, Phys. Rev. , 681 (1946) C. Kittel, Introduction to Solid State Physics (Wiley, New York,1996) G. Shvets, S. Trendafilov, J. B. Pendry, and A. Sarychev, Phys.rev. Lett. , 53903 (2007) F. Goos and H. H¨anchen, Annalen der Physik , 333 (1947) A. P. Hibbins, B. R. Evans, and J. R. Sambles, Science , 670 (2005) F. J. Garcia-Vidal, L. Martin-Moreno, and J. B. Pendry, Journal ofOptics A: Pure and Applied Optics , S97 (2005) S. A. Maier, M. L. Brongersma, and H. A. Atwater, App. Phys.Lett. , 16 (2001) B. D. F. Casse, W. T. Lu, Y. J. Huang, E. Gultepe, L. Menon, andS. Sridhar, App. Phys. Lett. , 023114 (2010) F. Lemoult, G. Lerosey, J. de Rosny, and M. Fink, Phys. Rev. Lett. , 173902 (2009) C. Lesueur, Rayonnement acoustique des structures (Editions Ey-rolles, 1988) F. Fahy and P. Gardonio, Sound and structural vibration: radia-tion, transmission and response (Academic Press, 2007) E. Yablonovitch, Phys. Rev. Lett. , 2059 (May 1987) P. de Vries and A. Lagendijk, Phys. Rev. Lett. , 1381 (Aug1998) M. Fleischhauer, Phys. Rev. A , 2534 (Sep 1999) S. Scheel, L. Kn¨oll, D.-G. Welsch, and S. M. Barnett, Phys. Rev.A , 1590 (Aug 1999) P. R. Berman and P. W. Milonni, Phys. Rev. Lett. , 053601 (Feb2004) M. Silveirinha, IEEE Transactions on Antennas and Propagation , 1766 (2006) E. Hendry, A. P. Hibbins, and J. R. Sambles, Phys. Rev. B78