OOrigin of folded bands in metamaterial crystals
Peter Markoˇs and Richard Hlubina
Department of Experimental Physics, Comenius University, Mlynsk´a Dolina F2, 842 48 Bratislava, Slovakia
Recently it has been found numerically that the spectra of metamaterial crystals may contain pairsof bands which disappear inside the Brillouin zone. We observe that the wave equations for suchsystems are essentially non-Hermitian, but PT -symmetric. We show that the real-frequency spectracorrespond to PT -symmetric solutions of the wave quation. At those momenta in the Brillouinzone where apparently no solutions exist, there appear pairs of complex-frequency solutions withspontaneously broken PT symmetry. PACS numbers: 02.60.Lj, 11.30.Er, 42.70.Qs, 78.67.Pt
One of the basic characteristics of waves propagatingin material media is their frequency spectrum. In peri-odic systems, for instance in photonic crystals, the fre-quency ω is a function of the wave vector (cid:126)q , ω = ω s ( (cid:126)q ),where (cid:126)q is restricted to an elementary tile of the recipro-cal space, the so-called first Brillouin zone. The discreteindex s numerates distinct branches of the dispersion,which correspond to different distributions of the wavefield within the unit cell of the periodic system. Sincein macroscopic systems the wave vector (cid:126)q changes quasi-continuously, each branch s leads in general to a finite in-terval of allowed frequencies, the so-called bands, whichmay be divided by band gaps in between them [1].In most systems studied so far, either in the solid-stateor photonic context, for each of the branches the function ω = ω s ( (cid:126)q ) stretches throughout the whole Brillouin zone.This is a simple consequence of Hermiticity. In fact, theplane wave e i(cid:126)q · (cid:126)x may experience a Bragg scattering toany of the plane waves of the form e i ( (cid:126)q + (cid:126)K ) · (cid:126)x , where (cid:126)K is a reciprocal lattice vector. In the basis of such states,the Schr¨odinger or wave equation takes the form of aneigenvalue problem H (cid:126)K (cid:126)K (cid:48) ( (cid:126)q ) c (cid:126)K (cid:48) = λ ( (cid:126)q ) c (cid:126)K . For a fixedcut-off we are then dealing with an N × N matrix which,if it is Hermitian, is guaranteed to have N real eigen-values, independently of the value of (cid:126)q . Smooth changesof H (cid:126)K (cid:126)K (cid:48) ( (cid:126)q ) lead then to smooth changes of λ ( (cid:126)q ), result-ing in bands which cannot disappear inside the Brillouinzone. In other words, the number of eigenfrequenciescannot be reduced in a certain interval of wave vectors.However, in numerical simulations it has recently beenfound that in certain systems the bands may disappearinside the Brillouin zone, forming the so-called foldedbands [2]. In particular, such behavior has been observedin photonic crystals in the form of a square array of meta-material cylinders immersed in vacuum. As an example,in Fig. 1 we show the two lowest-frequency bands forsuch a metamaterial photonic crystal calculated numeri-cally from transmission spectra [3]. Folded bands appearwhen the radius of cylinders R increases above the crit-ical value R c ≈ . a , where a is the spatial periodof the crystal. This surprising result indicates that thewave equation for electromagnetic field in a metamaterial photonic crystal is non-Hermitian [4].Actually, it is well known that if both, permittivity ε and permeability µ , are non-constant functions of thespatial coordinate, then the wave equation for the mag-netic field (cid:126)H reads M (cid:126)H ≡ µ − rot (cid:104) ε − rot (cid:126)H (cid:105) = ω (cid:126)H, (1)where the operator M is non-Hermitian even in ordinaryphotonic crystals made from dissipationless components[1, 5]. Note that we set the speed of light c = 1.So how can it be that folded bands have not been ob-served in ordinary photonic crystals? The reason is thatthe non-Hermitian character of Eq. (1) is often not es-sential, since it can be avoided by a reformulation of theproblem. For instance, if the permeability µ is real andpositive definite, one can redefine the magnetic field by (cid:126)H = (cid:126)h/ √ µ , thereby transforming Eq. (1) to the form O (cid:126)h ≡ µ − / rot (cid:104) ε − rot (cid:16) µ − / (cid:126)h (cid:17)(cid:105) = ω (cid:126)h (2)with an explicitly Hermitian operator O [1]. a q / a / R = 0.27 aR = 0.28 a π ω π aE z H x FIG. 1: (Color online) Dispersion relation ω = ω ( (cid:126)q ) in theΓ X direction for a two-dimensional photonic crystal made ofmetamaterial cylinders with real and dispersionless permit-tivity ε = − . µ = −
5, see inset. Electricfield is taken to be parallel to the cylinders. The two sets ofcurves correspond to cylinder radii slightly above and belowthe critical radius R c ≈ . a [3]. a r X i v : . [ phy s i c s . op ti c s ] O c t ω /2π -2-1012 f ( ) ql / ω / π ql / -0.200.20.40.60.8 ω / π π RealImag π l l l (a) (b) (c) ω FIG. 2: (Color online) 1D model with ε a = 1, µ a = 1, (cid:96) a =1, and ε b = − µ b = − (cid:96) b = 0 .
41. (a) Function f ( ω ).(b) Frequency spectrum. (c) Real and imaginary parts offrequency in the lowest folded band. The observation of folded bands in metamaterial pho-tonic crystals therefore suggests that their non-Hermitiancharacter should be essential , i.e. not avoidable by anyreformulation. In particular, it will be shown later that,in presence of interfaces between ordinary dielectric re-gions where √ µ is purely real and metamaterial regionswith purely imaginary √ µ , the operator O remains non-Hermitian [6].The goal of the present paper is to demonstrate thatthe appearance of folded bands in metamaterial pho-tonic crystals is a direct consequence of their essentialnon-Hermiticity. To this end, we will start by study-ing the simplest possible crystal structure, namely a one-dimensional (1D) periodic stack of right- and left-handedmaterials. Several anomalous features of electromagneticwave propagation have already been observed in thismodel [7–10]. Here we observe that the 1D model ex-hibits the so-called PT symmetry [11] and, making useof this recently developed concept, we will explain thepresence of folded bands in the spectrum of this model.Similar reasoning will be later applied to two-dimensional(2D) metamaterial crystals studied in [2, 3].We assume that the 1D stack consists of materials a and b characterized by ε i , µ i , refractive indices n i = √ ε i µ i , impedances Z i = (cid:112) µ i /ε i , and thicknesses (cid:96) i ,where i = a, b . All material parameters are assumed tobe real and frequency-independent. The frequency spec-trum of transverse electromagnetic waves, which propa-gate perpendicularly to the slabs with wave vector q , canbe determined from the implicit equation [12] f ( ω ) ≡
12 ( A +1) cos ωτ + −
12 ( A −
1) cos ωτ − = cos q(cid:96), (3)where (cid:96) = (cid:96) a + (cid:96) b is the length of the unit cell, τ ± = τ a ± τ b with τ i = n i (cid:96) i , and A = ( Z a /Z b + Z b /Z a ) / > a and b .In ordinary photonic crystals with ε i > µ i > τ + > | τ − | ≥
0. Therefore the larger-amplitude firstterm of f ( ω ) oscillates faster than the smaller-amplitudesecond term. Let us denote the positions of local extrema of the function f ( ω ) as ω ∗ . One can check easily [13]that | f ( ω ∗ ) | ≥ a is an ordinary di-electric with ε a > µ a >
0, whereas the slab b ismade from a metamaterial with ε b < µ b <
0, and n b <
0. In this case τ − > | τ + | ≥ f ( ω ) which oscillates fasterthan the larger-amplitude first term. As shown explic-itly in Fig. 2(a), then the values of f ( ω ∗ ) may lie withinthe interval ( − , A, β = τ − /τ + ) of Eq. (3).Since Eq. (3) does not shed light on the mathematicalstructure of the 1D problem, we shall restate its basicequations. The wave equation Eq. (1) reads piecewise − n − i H (cid:48)(cid:48) = ω H. (4)If we take the center of the slab a as the origin, then theboundary conditions at the interfaces ξ = ± (cid:96) a / H ( ξ − ) = H ( ξ + ) , E ( ξ − ) = E ( ξ + ) , (5)where ξ − and ξ + are infinitesimally shifted from ξ to theleft and right, respectively, and E ( ξ ± ) = H (cid:48) ( ξ ± ) /ε ( ξ ± ).Moreover, from the Bloch theorem follows an additionalboundary condition H ( (cid:96)/
2) = e iq(cid:96) H ( − (cid:96)/ . (6)The boundary-value problem which we have to solve isdefined by Eqs. (4,5,6).Let us prove that the operator O for the 1D prob-lem defined by Eqs. (4,5,6) is not Hermitian. To thisend, let us define magnetic and “electric” fields corre-sponding to h ( x ), H h ( x ) = h ( x ) / (cid:112) µ ( x ) and E h ( x ) = H (cid:48) h ( x ) /ε ( x ). If O was Hermitian, then the quantity D = (cid:82) (cid:96)/ − (cid:96)/ dx (cid:2) g ∗ ( x ) O h ( x ) − ( O g ( x )) ∗ h ( x ) (cid:3) should bezero for any pair of functions h ( x ) and g ( x ), both ofwhich generate fields H h , E h , and H g , E g satisfying theboundary conditions Eqs. (5,6). Integrating per partsand taking proper care of the boundary terms, we findthat D = ∆( − (cid:96) a /
2) + ∆( (cid:96) a / ξ ) = H h ( ξ ) [ E g ∗ ( ξ − ) − E g ∗ ( ξ + )]+ E h ( ξ ) [ H g ∗ ( ξ + ) − H g ∗ ( ξ − )]are the interface contributions. Note that in a right-handed medium [ H g ( x )] ∗ = H g ∗ ( x ) and [ E g ( x )] ∗ = E g ∗ ( x ), whereas in a left-handed medium [ H g ( x )] ∗ = ll ab | H | ( a r b it r a r y un it s ) -5-4-3-2-101 U + q l / -1-0.8-0.6-0.4-0.20 S spatiallydependent S complex ω π (b)(c) FIG. 3: (Color online) Observables for the 1D model withthe same parameters as in Fig. 2. (a) Spatial distribution of | H ( x ) | . White and shaded regions correspond to layers a and b , respectively. The vertical dashed lines denote the centers ofthe layers. Field distributions for the lowest folded band areshown for various values of cos q(cid:96) , from bottom to top: 0.650,0.700, 0.710, 0.719 (critical value), 0.730 and 0.800. Data areverticaly shifted for clarity. In PT -symmetric solutions thefield | H ( x ) | is even with respect to the center of any layer. PT symmetry-breaking solutions lack this symmetry [18]. (b) Theintegral quantity U + as a function of q . Note that U + = 0for PT symmetry-breaking solutions. (c) Real part of thePoynting vector, S , as a function of q . In all panels we usethe normalization | H ( (cid:96) a / | = 1. − H g ∗ ( x ) and [ E g ( x )] ∗ = − E g ∗ ( x ). It follows that inan ordinary photonic crystal the boundary conditionsEq. (5) imply that E g ∗ ( ξ − ) = E g ∗ ( ξ + ) and H g ∗ ( ξ + ) = H g ∗ ( ξ − ). Therefore D = 0, as was to be expected.However, in a metamaterial photonic crystal E g ∗ ( ξ − ) = − E g ∗ ( ξ + ) and H g ∗ ( ξ + ) = − H g ∗ ( ξ − ). Thus D (cid:54) = 0 andhence the operator O is non-Hermitian.Therefore the next question we should ask is: if theboundary value problem Eqs. (4,5,6) is non-Hermitian,how can its spectrum be real at all?In order to answer this question, let us define opera-tors of “parity” P and “time-reversal” T by ( P H )( x ) = H ( − x ) and ( T H )( x ) = H ∗ ( x ), respectively. We findthat Eqs. (4,5) are invariant under the action of both P and T , while Eq. (6) is invariant only under the com-bined antilinear operator PT . Therefore the 1D modelis PT -symmetric [15], as also observed for a somewhatsimilar boundary-value problem studied in [16].Since P = 1 and T = 1, the eigenvalues of the an-tilinear operator PT are λ = e iϕ , where ϕ is real [11].Therefore a PT -symmetric solution should satisfy H ∗ ( − x ) = e iϕ H ( x ) . (7)If the PT symmetry is not broken, i.e. if the eigenstatesof the wave equation are simultaneously also eigenstatesof PT , the eigenvalues ω have to be real [11]. On theother hand, if the eigenstates of the wave equation spon-taneously break the PT symmetry, then the eigenfre-quencies corresponding to H ( x ) and H ( x ) = PT H ( x ) a / M ω π Γ X M Γ XM FIG. 4: Dispersion relations for the 2D model with cylinderradius R = 0 . a along special lines in the 2D Brillouin zone(see inset). Electric field is polarized along the cylinders. Notethat in the Γ M direction, the two lowest bands are alreadyfolded. In the XM and Γ M directions, a third band entersthe studied frequency range. become complex conjugate. This is explicitly shown inFig. 2(c), which has been obtained by solving Eq. (3) un-der the assumption of a complex frequency ω = ω + iω .Such pairs of eigenvalues have been observed previouslyin the 1D problem [9] and also in other PT -symmetricproblems [17]. Note that the complex-frequency solu-tions appear exactly for those values of momentum q ,where the band folding has been observed. Therefore, ifwe allow for complex-valued frequencies, the number ofeigenvalues does not vary with q .We have checked numerically that for real-frequencysolutions, the magnetic field does obey Eq. (7). On theother hand, according to Eq. (7), a PT -symmetric so-lution should satisfy | H ( − x ) | = | H ( x ) | . But Fig. 3(a)clearly shows that this criterion is not satisfied forcomplex-frequency solutions, thereby proving explicitlythat they break the PT symmetry.Finally, we will turn back to the case of the squarearray of metamaterial cylinders with lattice constant a .Let us study electromagnetic waves with the electric fieldparallel with the cylinders and the z axis. Then thewave equation reads piecewise n − i (cid:52) E = ω E with a 2DLaplacian (cid:52) . The boundary conditions require the conti-nuity of E ( x, y ) and E (cid:48) n ( x, y ) /µ ( x, y ) across the cylindersurfaces, where E (cid:48) n denotes a normal derivative. More-over, a Bloch wave with wave vector (cid:126)q = ( q x , q y ) shouldsatisfy the boundary conditions E ( a/ , y ) = e iq x a E ( − a/ , y ) ,E ( x, a/
2) = e iq y a E ( x, − a/ . (8)Let us define the following 2D generalizations ofthe operators of “parity” P and “time-reversal” T :( P E )( x, y ) = E ( − x, − y ) ( T E )( x, y ) = E ∗ ( x, y ). Sincethe geometrical structure is P -invariant and since the ma-terial parameters ε i , µ i and n i are taken to be real, thewave equation and the boundary conditions inside theunit cell are invariant under both, parity P and time- qa / a / qa / π π π ω R = 0.1a R = 0.3a
FIG. 5: (Color online) Dispersion relations along Γ X for the2D model with cylinder radii different from Fig. 1. The shadedregions highlight the frequency gap (with complex momentum solutions inside) in the left panel, and the momentum gap(with complex frequency solutions inside) in the right panel.Note the space - time duality of the two regions. reversal T . However, similarly as in the 1D examplestudied before, the Bloch boundary condition Eq. (8) ispreserved only under the action of the combined opera-tor PT . Consequently, our 2D boundary value problemis again invariant under the action of the operator PT .Once the non-Hermiticity and PT symmetry of the 2Dproblem have been established, the reality of the spec-trum in the PT -symmetric sector of eigenstates follows.In analogy with the 1D case, we conjecture that in thoseregions of the Brillouin zone where the folded bands dis-appear, there exist pairs of complex-frequency solutionsin the spectrum. Unfortunately, since such states can notbe seen in a numerical transmission experiment [3], wecan not present a direct proof of their existence.In the 2D problem, the critical value of the cylinderradius R c , above which folded bands are formed, maydepend on the direction in momentum space. In fact,Fig. 4 shows that for instance for R = 0 . a , the twolowest bands are already folded in the Γ M direction.It turns out that one can find a simple criterion for theappearance of folded bands. In fact, from the Maxwellequations in the frequency domain, rot (cid:126)E = iω (cid:126)B androt (cid:126)H ∗ = iω ∗ (cid:126)D ∗ , after multiplication by (cid:126)H ∗ and (cid:126)E , re-spectively, follows the identity iω u − + ω u + + div (cid:126)S = 0 , (9)where u ± = (cid:126)E · (cid:126)D ∗ ± (cid:126)B · (cid:126)H ∗ and (cid:126)S = (cid:126)E × (cid:126)H ∗ . Notethat this identity holds also in dispersive and/or lossymedia and it represents a generalization of Poynting’stheorem for harmonic fields [19]. In our case of losslessmedia the quantities u ± are both real. Let us integrateEq. (9) over the area of an elementary cell of the photoniccrystal. The Bloch boundary conditions Eq. (8) implythat (cid:82) cell d x div (cid:126)S = 0. Therefore iω U − + ω U + = 0 musthold, where U ± = (cid:82) cell d x (cid:104) ε | (cid:126)E | ± µ | (cid:126)H | (cid:105) . Consideringthe real part of this condition for PT symmetry-breaking solutions with ω (cid:54) = 0, we find that U + = 0 must hold,as confirmed for the 1D problem in Fig. 3(b). Similarcriteria for the occurence of folded bands have been foundin [2] by different reasoning. Note that at the verge of PT symmetry breaking, the group velocity diverges and U + vanishes. Explicit calculation for the 1D problem showsthat the real part of (cid:126)S stays finite here, see Fig. 3(c).Before concluding let us observe that folded bands areclosely related to the frequency gap in periodic systems.In fact, it is well known that within the frequency gap,there exist only solutions with a complex wave-vector.On the other hand, folded bands imply the existence ofa momentum gap , inside which there appear only stateswith a complex frequency. Thus the momentum gap is aspace - time dual of the frequency gap, see Fig. 5.In conclusion, we have explained the physical origin offolded bands, or momentum gaps, in spectra of metama-terial crystals. Two ingredients are necessary for theirappearance: first, the boundary-value problem has to beessentially non-Hermitian, and second, the possible re-ality of its spectrum has to be guaranteed by additionalsymmetry, such as the PT symmetry in the metamaterialcase. In systems fulfilling these assumptions, momentumgaps should be commonplace.We thank V. Balek and M. Mojˇziˇs for helpful discus-sions. This work was supported by the Slovak Researchand Development Agency under the contract No. APVV-0108-11 and by the Agency VEGA under the contractsNo. 1/0372/13 and No. 1/0904/15. [1] K. Sakoda, Optical Properties of Photonic Crystals,Berlin, Heidelberg: Springer (2005).[2] P. Y. Chen et al , New J. Phys. et al , Phys.Rev. B et al. , Photonic Crystals: Moldingthe Flow of Light, 2nd ed. Princeton: Princeton Univer-sity Press (2008).[6] Since we are dealing with a metamaterial crystal, Her-mitization of the wave equation for (cid:126)E in terms of thesubstitution (cid:126)E = (cid:126)e/ √ ε is plagued by the same problem.[7] I. S. Nefedov and S. A. Tretyakov, Phys. Rev. E ,036611 (2002).[8] Jensen Li et al. , Phys. Rev. Lett. , 083901 (2003).[9] Liang Wu, Sailing He, and Linfang Shen, Phys. Rev. B , 235103 (2003)[10] D. Bria et al. , Phys. Rev. E , 066613 (2004). [11] For a review, see C. M. Bender, Contemporary Physics , 277 (2005).[12] P. Yeh, Optical Waves in Layered Materials, Hoboken,New Jersey: Wiley (2005).[13] To this end it suffices to note that for frequencies ω n = nπ/τ + , where n is an integer, we have | f ( ω n ) | ≥
1. More-over, f ( ω n ) exhibits even-odd oscillations with n . Butsince the second term in f ( ω ) oscillates with a longer pe-riod, in between ω n and ω n +1 there will be at most oneextreme of f ( ω ), and therefore | f ( ω ∗ ) | ≥ | f ( ω n ) | ≥ n gaps [8] are associated witha special case of folded bands. In fact, the zero-¯ n condi-tion is equivalent to τ + = 0 and in this case Eq. (3) hassolutions only for q = 0 and ωτ − = 2 mπ , where m is aninteger.[15] In previous work, the concept of PT symmetry has beenexploited in optics from a quite different perspective,making use of the formal analogy between quantum me-chanics and the paraxial approximation. The PT symme-try was achieved by judicious spatial modulation of loss and gain, see e.g. C. E. R¨uter et al. , Nat. Phys. , 192(2010). The closest in spirit to ours is the paper L. Ge andH. E. T¨ureci, Phys. Rev. A , 053810 (2013), which alsostudies metamaterial crystals (with different symmetry)using the concept of PT symmetry. That paper is limitedonly to 1D systems, however, and, more importantly, itdoes not study the eigenvalue problem.[16] D. Krejˇciˇr´ık, H. B´ıla, and M. Znojil, J. Phys. A: Math.Gen. , 10143 (2006).[17] For a recent example, see A. Mostafazadeh, J. Phys. A:Math. Theor. , 375302 (2011).[18] In the PT symmetry-breaking sector, we plot | H ( x ) | for only one solution for every value of q . The minimaof | H ( x ) | shift to the left from the center of layer a .The other solution (not shown), which corresponds to acomplex-conjugate frequency, exhibits an opposite shift.[19] J. D. Jackson, Classical Electrodynamics, 3 rdrd