Beam oscillations and curling in chirped periodic structures with metamaterials

We study the propagation of electromagnetic waves in one-dimensional chirped periodic structures composed of alternating layers of negative-index (or left-handed) metamaterial and conventional dielectric, under the condition of the zero average refractive index. We consider the case when the periodic structure has a linear chirp, and the chirp is introduced by varying linearly the thickness of the layers across the structure. We apply an asymptotic analytical method for the analysis of the Bloch oscillations and find that, in a sharp contrast to ordinary periodic dielectric structures, the energy flow in multi-layered stacks with metamaterials may have the opposite direction at the band edges, thus providing novel possibilities for the beam steering in the transmission band.


I. INTRODUCTION
Materials with negative index of refraction, also known as left-handed metamaterials, are attracting a great scientific interest nowadays. This novel type of materials was predicted theoretically back to 1967 [1] but only 30 years later this theoretical curiosity was followed by the first experimental demonstrations based on the fabrication of composite structures consisting of split-ring resonators and metallic wires [2,3]. The basic properties of the left-handed metamaterials predicted theoretically have been confirmed by experiment, and this led to a rapid progress in this field and many other interesting discoveries. In particular, many new physical phenomena based on the concept of metamaterials and negative index of refraction have been predicted, including the fundamental concept of perfect lens [4] and, more recently, electromagnetic cloaking [5,6,7].
One of the interesting directions in the analysis of the unusual properties of this novel composite metamaterials is the study of periodic structures (or photonic crystals) composed of metamaterials. Photonic crystals, being widely used for both light control and beam manipulation, demonstrate a variety of novel physical effects, and they may offer many new possibilities being combined with metamaterials. In particular, a series of recent studies revealed the existence of a novel zero-index bangaps [8,9,10,11] associated with one-dimensional photonic crystals composed of alternating layers of negativeindex and conventional dielectrics for which the averaged refractive index vanishes. Such novel periodic structures demonstrate many intriguing properties, including substantial suppression of the Anderson localization and long-wavelength resonances [12].
In this paper, we study the propagation of electromagnetic weaves in one-dimensional periodic structures composed of alternating layers of left-handed metamaterial and conventional dielectric. It is well established [13,14] that introducing a chirp into a one-dimensional periodic structure leads to the generation of optical Bloch oscillations which can be also predicted in the presence of metamaterials [15]. In addition to our numerical analysis of the Bloch oscillations studied earlier [15], in this pa-per we demonstrate a number of novel, unique features of the beam propagation in this kind of chirped periodic composite structures. In particular, in contrast to Ref. [15], here we consider the chirped structures with a large number of layers and slowly varying period. Applying the approximation based on the geometric optics, we reveal a variety of novel effects that could be observed in such structures, including the reversal of the energy flow at the band edges and the fascinating effect of the beam curling. This novel effect of the beam curling can be useful for the beam steering in the transmission band. We also study the eigenvalue problem and the corresponding eigenmodes and determine both the field structure and Poynting vector distribution associated with the beam propagation.
The paper is organized as follows. In Sec. II we introduce our structure as a one-dimensional stack of layers composed of alternating layers of metamaterial and dielectrics. We discuss the way of presenting eigenmodes of the structure with a slowly varying period and discuss the initial conditions. In Sec. III we implement the Bloch theorem and derive the dispersion relation for the Bloch wave vector. Then, we analyze the bandgap diagram for this structure and discuss the conditions for the observation of a novel type of the beam dynamics associated with the Bloch oscillations. Section IV is devoted to the geometric optics approximation. Here we demonstrate a possibility to apply the geometric optics for the analysis of layered structures with slowly changing parameters. We derive the equations of motion for the paraxial beams, and also study the period of oscillations pointing out novel properties of the beam trajectories connected with a change of the sign of the spatial group velocity in the structure due to the presence of a metamaterial. Also, we analyze different regimes of the beam propagation and compare them with the well-known case corresponding to the similar structures composed of conventional dielectric layers. In Sec. V, we discuss the eigenvalue problem and point out the difficulties and failure of the commonly used numerical techniques. In addition, we describe an alternative method to find the eigenmodes of the structures with a large number of periods. Based on this method, we study the field distribution in the structure and determine the direction of the Poynting vector in the beam. We observe a close correspondence with predictions of the geometric optics. Finally, Sec. VI concludes the paper.

II. STRUCTURE GEOMETRY AND EIGENMODE FORMALISM
We study a one-dimensional chirped periodic structure shown schematically in Fig. 1, where the slabs of a negative-index metamaterial with the width d l m are separated by the layers of the conventional dielectric with the width d r m. We describe the variation of the refractive index in the m-th pair of layers as follows, where n l and n r are the refractive indices of metamaterial and dielectric layers, respectively, Λ m is the width of the m-th unit cell. In contrast to our previous studies [15], here we analyze the structures with a large number of unit cells and slowly varying period. We consider TE-polarized waves with the electric field described by one component E = (E x , 0, 0), and the waves propagating in the (y, z) plane. In this case the problem is described by the Helmholtz equation, where ∆ 2 is the two-dimensional Laplacian. The field is assumed to be monochromatic with the frequency ω, and the coordinates are expressed in the units of c/ω, c is the speed of light. For the structures with slowly varying periods, the eigenmodes of the problem (2) can be introduced by applying the Bloch theorem for periodic systems and the corresponding Bloch-wave formalism [16]. In this approach, the electric field corresponding to the eigenmode with the propagation constant k y in the unit cell with the width Λ can be represented as, where k y and K b are normalized to ω/c, the amplitude A Λ varies slow across the structure, U Λ is the Bloch function, and K Λ b is the Bloch number of the periodic structure with the unit cell size Λ, which we analyze systematically in Sec. III.
We study propagation of the Gaussian beams with the electric field at y = 0 defined as where E(z) is the eigenmode of the Eq. (3) with the propagation constant k y = k 0 y , which will determine propagation direction; h is the beam width, which we consider to be significantly large, i.e. h ≫ < Λ >. the structure can be represented as a superposition of eigenmodes. As long as the Gaussian beam is wide with respect to the unit cell, the spectrum of the contributing eigenmodes is narrow, and it is centered near the point k 0 y [16]. Such beam can be treated within the paraxial approximation. The beam propagation is mainly described by the modes with k y ≃ k 0 y and the tangent of the beam angle is defined by (∂k y /∂K b ) |ky=k 0 y ,Λ=Λ0 , where Λ 0 corresponds to the beam launching point. An additional beam tilt can be introduced in more general case of complex ψ(z).

III. PERIODIC STRUCTURE: BANDGAP PROPERTIES
As the first step of our analysis, we study the periodic structure with the period Λ without a chirp. In this case Eq. (2) has periodic coefficients, and we can apply the Bloch theorem [17] and present the electric field in the form, where K b is the Bloch number, and the field envelope U (z) is a periodic function with the period Λ, so that U (z + Λ) = U (z). In the m-th pair of layers, the Bloch waves are presented as, where the indices l and r correspond to the left-and righthanded slabs, respectively, and the amplitudes a l,r and b l,r are found from the boundary conditions at the interfaces separating the metamaterial and dielectric layers (see, e.g., Refs. [16,17,18]). The Bloch wavenumber K b is defined from the dispersion relation, where k zl,zr = ∓ n 2 l,r − k 2 y and k y is the propagation constant along the y axis.
According to Eq. (5), the waves can propagate in the structure when K b is real. In Fig. 2(middle), we plot the bandgap diagram of the layered structure on the plane (Λ, k y ). Here we assume that the dielectric layer is air, ε r = µ r = 1, and the dielectric layers are two times thicker than the metamaterial layers, i.e. d rm /d lm = 2, and this ratio is preserved across the structure. We choose the metamaterial with the parameters ε l = −5 and µ l = −0.8.
The bandgap properties of one-dimensional periodic structures with metamaterials have been studied in several papers [10,11], including the case of the zero-index averaged refractive index [9]. It was shown that for the structures with Λ n(z)dz = 0 the bandgap spectrum includes transmission resonances which shrink into infinitesimally thin lines into a complete bandgap when the widths of the slabs coincide [9]. For the normal incidence (i.e. when k y = 0), the transmission is observed only when n r d r = n l d l = πq, where q is integer.
We study the Bloch wave phase incursion K b Λ across the structure for different values of the propagation constant k y , see Fig. 2(left). For small values of the propagation constant (e.g., k y < 0.5), the low-order transmission resonances have zero phase incursion, in contrast to the case of conventional dielectric Bragg gratings. The maximum of the phase is accumulated in the middle of the transmission band, and its amplitude increases with the increase of the band order, reaching the value of π in the n-th order transmission resonance. A growth of the propagation constant leads to the band coupling (0.5 < ∼ k y < ∼ 1), not shown here. The phase incursion in a new coupled broadband region is approaching π. Further increase of the propagation constant leads to the narrowing of the transmission regions and to a rapid increase of the phase across the band to the value of π. On the bandgap diagram, we observe two turning points at k y ≃ 1.7 and k y ≃ 2, (marked points 1 and 2 in the inset of Fig. 3(middle)). At these points the modes with close k y have the same value of the Bloch wavenumber near the band edge. We trace the behavior of the band edges ζ − (k y ) and ζ + (k y ) with a change of the propagation constant k y . For k y ≃ 1.7, the left band edge ζ − changes the type of its monotonicity. Consequently, in this region the left (ζ − ) and right (ζ + ) band edges have different slopes, i.e. (∂ζ − /∂k y )(∂ζ + /∂k y ) < 0, as shown in the inset of Fig. 3. The latter, with consideration that phase is growing from 0 to π in the transmission region, leads to the crossing of the corresponding modes (see inset in Fig. 2). Finally, for k y ≃ 2, the right band edge ζ + also changes the sign of its slope, and (∂ζ − /∂k y )(∂ζ + /∂k y ) becomes positive. We note that for small values of k y , the first-order band edges also have different signs of their slopes but a new type of behavior cannot be observed because the phase incursion vanishes in the band for k y < 0.5. Also, in the conventional Bragg gratings both band boundaries have the same slope and thus the mode crossing does not occur.
Crossing of two modes means that the effective spatial velocity defined as ν = ∂k y /∂K b , has different signs at different sides of the crossing point, and it becomes singular at the point itself. As long as the velocity is related to the energy flow, we observe different directions of the energy flow in the y-direction at different edges of the band, as will be discussed in more details in Sec. IV and Sec. V. We plot the velocity for different values of the propagation constant k y , see Figs. 2(a-c). For k y outside the region, Fig. 2(a), the velocity does not change its sign in the band, thus meaning that energy flow has the same direction at any point of the band. Similar type of the beam evolution can be found in the stacks of conventional dielectric layers. Further increase of k y leads to more asymmetric velocity profile and, finally, at k y ≃ 1.7 [see Fig. 2(b)], the velocity becomes infinite at the left band edge. In Sec. IV below we demonstrate that this case corresponds to the energy flow vanishing in the y-direction. When the propagation constant k y is inside the region [see Fig. 2(c)], the velocity profile has the second-order discontinuity with the sign changed inside the band. Consequently, the energy flow in the y-direction also changes its sign within the transmission band. The singularity point in the spatial velocity profile moves from the left band edge to the right edge with the growth of the propagation constant from k y ≃ 1.7 to k y ≃ 2, (see turning points 1 and 2 in Fig. 3). Finally, for the propagation constants k y > 2, the velocity profile becomes continuous and negative, indicating that the waves are backward in the whole band.
Below, we consider the simplest case of a linear chirp of the structure, Λ i+1 = Λ i + δΛ with δΛ ≪ Λ, see Fig. 1. For δΛ → 0, i.e. an adiabatic change of the unit cell size across the structure, we implement the geometric optics approximation.

IV. GEOMETRIC OPTICS APPROXIMATION
We study the evolution of the Poynting vector by employing an analogy with a homogeneous medium having gradually changing refraction index [18,19]. We consider the time-averaged Poynting vector, Re where E x is given by Eqs. (4), e y and e z are the unit vectors. The time-averaged electromagnetic energy density is defined as Re(εEE * + µHH * ) = (10) Re In a homogeneous medium, the time-averaged Poynting vector is tangential to the beam trajectory. It changes gradually with continuously varying refractive index, thus defining a geometric ray. In our case, the Poynting vector changes its direction stepwise within each unit cell, since the refractive indices of the slabs comprising the unit cell have the opposite signs. To describe an average behavior of the energy flow on larger scale, we consider the Poynting vector and energy density averaged We introduce the average velocity of the energy flow [16,17]: As long as the width of the unit cell changes adiabatically across the structure, the average velocity of the energy flow v e changes quasi-continuously. Therefore, we can define the rays in the space (Λ, y) (we note here that due to the linear chirp of the structure, the spatial coordinate z m is a linear function of Λ) as trajectories along which the average energy < S > is directed. Thus, the ray equation is defined as follows: where r = e y y + e z Λ is the radius-vector describing the ray. In Ref. [16] it was shown that for a paraxial beam the average velocity of the energy flow coincides with the group velocity defined as where K = e y k y + e z K b is the wave vector, and ω is the angular frequency. For a fixed frequency ω of monochromatic beam we obtain the dispersion relation k y = f (K b ), (see Eq. 7), corresponding to the curve in the space (k y , K b ) described by the wave vector K, which is analogous to the equifrequency surface for anisotropic media [20]. Consequently, using Eqs. (11), (12), and (13) we obtain v g = v e and < S > are normal to the curve k y = f (K b ). The latter condition is described by the equation, where < S y > and < S z > are the y and z components of the Poynting vector < S >, ν is the spatial velocity of the ray. Relation (14) between the Poynting vector components and velocity ν shows that when ν → ∞, < S y >= 0, and the total energy flows along the z-direction. On the other hand, when ν = 0, the energy flows along the structure (< S z >= 0).
Using Eqs. (11) and (14), we rewrite Eq. (12) as follows: Now we can derive the ray trajectory equation in the space (Λ, y) where ∓ correspond to the forward and backward propagating waves, respectively. We note that the motion equation (16) can also be derived from the Fermat principle. We multiply Eq. (12) period of oscillations, Φ by K and using Eq. (13), we write the Fermat principle in the form [20]: where dr = e y dy + e z dΛ is the unit vector along the beam trajectory. The Fermat principle is the principle of least action with the coordinate y playing the role of time. The Lagrange function L = ±K b (dΛ/dy) + k y is defined accurate to the full y-derivative of Λ. As a result, the unambiguity in the Bloch wavenumber, K b = K (0) b + 2πm/Λ does not affect the ray behavior. The Hamilton equations for the rays can be derived in the form [19]: We remind now that k y remains constant across the structure, consequently ∂k y /∂Λ = 0, meaning that the problem is invariant to translation along the layers. Hence, we consider only the first equation, which can be also derived from general suggestions, assuming that k y refers to the "energy", being preserved in the structure and ∂k y /∂K b being the "group velocity" of the paraxial beam.
The first equation of the system (18) describes diffractionless motion of a beam with a narrow spectrum centered near the point k 0 y . The beam motion is equivalent to the motion of an effective particle in a one-dimensional potential W = [ν 2 (Λ)/2] between left (ζ − ) and right (ζ + ) band edges. As long as the spectrum of allowed energies near the point k 0 y is quasi-continuous for infinite structures with adiabatically changing width Λ, the motion remains periodic with the period of oscillations Φ found as, where ζ − and ζ + refer to the band edges, ν is defined by Eq. (14). In the gap regions, the field decays exponentially across the structure so that the motion is prohibited, and the spatial velocity vanishes at the band edges. Consequently, the expression under the integral has singular points at ζ ± , and the integration fails. To resolve this problem, we introduce a new variable ξ, and obtain The expression under the integral is finite within the integration boundaries, hence the new integral is converging. Using the integral (21) we calculate the dependence of the period on the propagation constant k y , see Fig. 3(left). We observe that with a growth of k y the oscillation period grows as well, this is explained by broadening of the band described in Sec. II. A further increase of k y leads to the beam narrowing and decrease of the oscillation period Φ. For 1.7 < k y < 2, i.e between the turning points, the period vanishes meaning that the corresponding mode does not transfer energy along the structure. A further increase of k y leads to negative values for the period of oscillations, which corresponds to the total energy flow in the direction opposite to the propagation constant k y .
Using Eq. (18) and Eq. (21), we calculate the beam trajectories. First of all, we verify our results for the well-known structures with conventional dielectric slabs, and the trajectory corresponding to this case is shown in Fig. 4(left). We assume that a dielectric layer with the width b has the parameters ε b = 10 and µ b = 1, and the layers are separated by vacuum of the width a, so that a/b = 2. We reveal a close correspondence between our results and the trajectories calculated by means of the geometric optics in Ref. [14]. The beam is reflected from both band boundaries, where the spatial velocity vanishes. The total energy flow along the y-direction is positive, and it preserves its sign during oscillations. For the structures with metamaterials, the beam with the spectrum centered near the propagation constants k y outside of the region between turning points (points 1 and 2 in Fig. 3) experiences the same reflection from the boundaries of the Brillouin zone with the total energy flow along the structure in the positive direction, for k y < 1.7, and in opposite (negative) direction, for k y > 2.
Next, we study the beam propagation for the values of k y inside the region between the turning points. The trajectory calculated near the first turning point (see Fig. 4) has a peculiarity near the left band edge. Corresponding period of oscillations is shown in Fig. 3(left). Near the left band edge the group velocity becomes infinite, thus energy flow along the layers vanishes, and the energy flow across the structure changes its sign. Hence a vortex structure is formed at this point. The observed spike is formed by narrowing of the trajectory near the left band edge with increase of propagation constant until it reaches the first turning point. Further increase of propagation constant, as was described in Sec. II, corresponds to a shift of the singular point of the group velocity further into the band region. The trajectories corresponding to this case are shown in Figs. 4(a-e). Since the group velocity changes its sign before and after the singular point, the trajectory crosses itself and a beam curl is formed. For k 0 y = 1.85 we observe a curl near the left band edge of the structure, see Fig. 4(b). Further increase of the propagation constant leads to the increase of the curl's size, and decrease of the period, see Fig. 3. Finally, at k 0 y ≃ 1.895 the curl reaches the other transmission band edge forming a practically closed trajectory with almost zero period of oscillations, see Fig. 4(c). Note that near the left band edge the energy flow in the y-direction is negative but remains positive near the right band edge. Further increase of the wavenumber results in the opposite process with the curl forming near the right band edge, see Fig. 4(d). For the propagation constants k y corresponding to the second turning point the curl vanishes, with a spike forming near the right band edge, and the total energy flow along the structure is in the opposite direction to the propagation constant. For k 0 y > 2, the spike vanishes, and smooth reflection from the right band edge is formed. The trajectory again becomes continuous and without any ambiguities.

V. NUMERICAL SIMULATIONS
For the finite structures with non-adiabatic change of the unit-cell size, the approximation of geometric optics is not valid, however we can still apply the Bloch-wave formalism for describing the beam propagation, including the beam diffraction and formation of interference patterns. To find the field distribution in the structure, we calculate the eigenmodes E i (z) and eigenvalues k i y of a finite stack of layers. There are several known approaches to solve this problem. The first approach is based on the representation of the field in each slab in the form of the counter-propagating waves. Boundary conditions between the layers define the ratio of the corresponding wave amplitudes in the slabs. These amplitudes can be represented in terms of either transfer matrix or scattering matrix formalisms [17,21]. The transfer matrix approach can be used for calculating the eigenmodes and eigenvalues for the structures with a small number of periods, and it is not suitable for calculating evanescent modes in long enough structures. The scattering matrix approach provides a better convergence, but it also cannot be applied for calculating evanescent modes in the structures with more than approximately 100 slabs. Another approach is based on the discretization of the wave equation across the structure, and we have used it previously to find the eigenvalues of the discrete problem [15]. However, this approach is computationally intensive, and it is not practical for calculating the transmission of very long structures.
Here we implement another approach based on the Bloch wave formalism. In Sec. II we have already demonstrated that for the slowly varying period of the structure (i.e., δΛ ≪ Λ) the eigenmodes can be described by Eq. (3) with the amplitudes A Λ and wavenumbers k y . We assume that the structure with a linearly changing period is placed between two semi-infinite periodic structures, with the period coinciding with that of the adjacent layers of the structure, see Fig. 1. In each unit cell the field can be presented as a superposition of the forward and backward propagating Bloch waves, where A Λ and B Λ are the slowly varying amplitudes, U ± and K b are defined by Eqs. (6) and (7), respectively.
In the left semi-infinite structure, the fields decay exponentially, and we choose the amplitudes in the first layer as A Λ l = 1 and B Λ l = 0. In the right semi-infinite structure the forward wave vanishes, and this means that A Λr = 0 and B Λr = 1. We calculate the amplitudes in the adjacent unit cells using the boundary conditions. Repeating this procedure we find the amplitudes in the whole structure (for more details, see Appendix A). We note that K b is real in the band, hence the Bloch waves are propagating, and the modes do not grow exponentially. In a small number of layers adjacent to the band, K b is imaginary and, consequently, the amplitudes grow exponentially. We choose relatively small number of unit cells with the parameters corresponding to the bandgap, so that the amplitudes A Λ and B Λ remain reasonable across the structure. We note that the absolute value of the ratio of the amplitudes of the counter-propagating waves in the band is unity, so that the energy is localized in the transverse direction inside the structure. We trace the dynamics of the Bloch wave amplitudes from left and right boundaries separately and compare the amplitudes in the middle of the band region, see Fig. 5. For the eigenmodes with the eigenvalues k y , the amplitudes of the waves at the left and right edges of the structures become linearly dependent, i.e., Using the procedure described above, we find the propagation constants k y corresponding to the eigenvalues and eigenmodes of the problem. Then, the initial field distribution is represented as a superposition of the eigenmodes, and the full structure of the fields can be retrieved. We calculate the eigenmodes for the structure consisting of 2100 unit-cells with the width changing linearly from the value Λ l = 4.23 to the value Λ r = 4.44, and δΛ ≃ 2×10 −4 . We note that < Λ >= (Λ r +Λ l )/2 ≃ 2/3λ, where λ is a free space wavelength. The width of the whole structure is approximately 1500λ, which is, e.g., about few millimeters in infrared regime. We have chosen the initial field distribution in the form of Gaussian beam, see Eq. (4), launched near the right band boundary of the structure, see k 0 y = 1.838. We decompose the initial field distribution into the superposition of eigenmodes of the structure by the least-squares method. The spectrum of the excited eigenmodes is shown in Fig. 3(right). The spectrum is narrow, and it is centered near k y = 1.838. The spectrum is to some accuracy equidistant, meaning that the field restores its shape after the distance Φ = 2π/δk y . We plot the field amplitude averaged over the unit cell. The beam curling predicted in the framework of the geometric optics is clearly observed. At the beam self-crossing point we observe an interference pattern created by the forward and backward propagating waves. The field distribution is restored almost completely after the first period of oscillations. Now we calculate the distribution and direction of the Poynting vector. In Fig. 7 we plot the Poynting vector averaged across the unit cells. It is clearly seen that near the left boundary the direction of the Poynting vector along y-axis is opposite to the direction near the right boundary. At the self-crossing point corresponding to singularity of the group velocity, we observe that the energy flow vanishes in the y-direction.
We also calculate the average y-component of the Pointing vector, defined as < S y >, on the unit cell and show its dynamics across the structure for the mode with the propagation constant k y = 1.838, see Fig. 5. It is clearly seen that at the left boundary of the band the energy flow direction is opposite to that near the right boundary, and it vanishes in the center of the structure. This dynamics is consistent with the behavior predicted within geometric optics approximation in Secs. II and IV. We note that the total energy flow, i.e. the integrated flow across the structure, remains positive for this mode and for all modes contributing to the beam spectrum.

VI. CONCLUSIONS
We have presented a systematic analysis of the propagation of electromagnetic waves in chirped periodic structures composed of two kinds of alternating layers, the layers of negative-index metamaterial and the layers of conventional dielectrics, under the condition of the zero averaged refractive index. We have considered the chirp in the structure parameters introduced by varying the thickness of all layers linearly across the structure, and we have applied the methods of the geometric optics for analyzing the beam propagation and Bloch oscillations in such infinite composite structures. For the adiabatically changing period of structure, we have predicted the beam self-crossing in the bands, and have found that the energy flow in such multi-layer structures with metamaterials may have the opposite direction at the band edges, in a sharp contrast to the similar structures composed of conventional dielectrics. This novel effect of the beam curling can be useful for the beam steering in the transmission band.