**Z**-transform is developed to verify the analytical results. The approach can be applied to arbitrary periodic refractive index profiles using the Fourier series method.

## 1.

## Introduction

Electromagnetic metamaterials (MM) are artificial composites with electromagnetic properties not readily found in nature. A special class of MMs is the negative refractive index metamaterials (NRM), artificial structures with negative phase velocity.^{1} A number of practical implementations of optical MM have been reported.^{2}^{,}^{3}

NRM are typically produced using arrays of subwavelength “particles” with negative effective relative permittivity and permeability as their structural units. The first proposed NRM particles were split-ring resonators and nanowires, furnishing negative permeability and permittivity of their composites.^{3} They are well understood and extensively used in the microwave domain. However, other particles such as complementary split-ring resonators,^{4} cut-wire pairs/plate pairs,^{5} and double ﬁshnets^{6}7.^{–}^{8} are also investigated. The first NRM were experimentally confirmed in 2001,^{9} and recently the experimental fishnet-type NRM in the visible range of frequencies have been fabricated and investigated.^{10}

The properties of NRM, such as the negative index of refraction (and negative phase velocity), inverse Doppler effect, radiation tension instead of pressure, etc.^{11}^{,}^{12} resulted in a number of proposed applications. Among those we mention superlenses and hyperlenses that enable imaging far below the diffraction limit,^{13}^{,}^{14} resonant cavities, and waveguides with geometrical dimension orders of magnitude smaller than the operating wavelength^{15} as well as invisibility cloaks and generally transformation optics.^{16}

Most studies consider structures with constant effective permittivity and permeability within the NRM part and abrupt interfaces with the surrounding regular materials [positive refractive index media (PRM)]. However, there is both theoretical and practical interest in NRM with spatially varying effective permittivity and permeability within the NRM structure and with gradual transition from the PRM to NRM and vice versa. Graded refractive index is interesting for transformation optics including hyperlenses^{17} and invisibility cloaks.^{18} Various other proposed applications of gradient refractive index (GRIN) metamaterials include beam shaping and directing, enhancement of nonlinear effects,^{19} superlenses,^{20} etc.

The first paper dedicated to GRIN NRM was published in 2005.^{21} Analytical approaches to graded index metamaterial structures are of special interest, since they ensure fast, simple, and direct route to the determination of the field distribution and the calculation of the scattering parameters within such materials.^{22}23.24.25.26.27.^{–}^{28}

The present study is a generalization of our previous work^{27}^{,}^{28} in a sense that it allows different loss factors in PRM and NRM segments. An assumption of uniform loss factors throughout the structure has been made in Refs. 2728.–29. The possibility of choosing arbitrary loss factors in PRM and NRM, independent of each other, is essential for a realistic description of composites involving NRM as building blocks, since it is empirically well known that the losses in the NRM part are significantly higher than those in the PRM. Thus, in this paper, we present an exact analytical solution of Helmholtz equations for the propagation of electromagnetic waves through a periodic gradient-index PRM–NRM composite with most general loss factors in the two materials and with sinusoidal periodicity for the case of constant impedance throughout the structure. A comparison of the obtained analytical solution to the results of numerical simulation using a $Z$-transform based model is given.

## 2.

## Problem Formulation

The geometry of the present problem is shown in Fig. 1. The electric field points to the $y$-direction and has the form $\overrightarrow{E}(\overrightarrow{r})=E(x){\overrightarrow{e}}_{y}$, while the magnetic field points to the $z$-direction and has the form $H(\overrightarrow{r})=H(x){\overrightarrow{e}}_{z}$. The wave propagates along the $x$-axis. The spatial variation of the refractive index along the $x$-axis is described by a cosine function. An implementation of a GRIN metamaterial obtained following an approach similar to that found in Ref. 29 is shown in Fig. 2. For this purpose, a two-dimensional array of single split rings is deposited onto a dielectric substrate. The example is given for illustration purposes only, as the structures with varying dimensions may be any of metamaterial unit structures or “particles” (“atoms”). Also, periodicity does not have to be sinusoidal (and actually may be represented by any graded structure, as long as a Fourier series representation is valid) and the gradient itself may be along one, two, or all three axes.

The NRM structures are of importance for transformation optics,^{30} the most well-known example being the invisibility cloaks. For instance, optical carpet cloaks^{31} were reported with effective index gradient obtained by drilling hole arrays with varying geometry.^{32} Other important applications include optical and generally electromagnetic concentrators based on metamaterials, beam shapers, and beam steering devices, as well as different kinds of metamaterial lenses including the hyperlenses^{17} for the transformation of near field into the far field. Finally, an important application is gradient index circuitry utilizing metamaterial waveguides.^{33}

We write the Helmholtz equations for $E(x)$ and $H(x)$^{28}

## (1)

$$\frac{{\mathrm{d}}^{2}E}{\mathrm{d}{x}^{2}}-\frac{1}{\mu}\frac{\mathrm{d}\mu}{\mathrm{d}x}\frac{\mathrm{d}E}{\mathrm{d}x}+{\omega}^{2}\mu \epsilon E(x)=0,\phantom{\rule[-0.0ex]{2em}{0.0ex}}\frac{{\mathrm{d}}^{2}H}{\mathrm{d}{x}^{2}}-\frac{1}{\epsilon}\frac{\mathrm{d}\epsilon}{\mathrm{d}x}\frac{\mathrm{d}H}{\mathrm{d}x}+{\omega}^{2}\mu \epsilon H(x)=0,$$## 3.

## Analytical Solutions of the Field Equations

Let us now consider an infinite and inhomogeneous periodic structure, where the real parts of the effective dielectric permittivity and magnetic permeability vary as a cosine function from positive values (right-handed material) to negative ones (left-handed material) and back again. The thickness $a$ of the positive part is equal to that of the negative part. For the sake of simplicity, we assume an impedance-matched case where real parts of effective permittivities and permeabilities of the two materials at a given frequency have opposite signs and equal absolute values. Thus, to determine their values, we use the functions

## (2)

$$\mu (\omega ,x)={\mu}_{0}{\mu}_{R}(\omega )\mathrm{cos}\left(\frac{\pi x}{a}\right)-i{\mu}_{0}[\frac{{\mu}_{\mathrm{IR}}+{\mu}_{\mathrm{IL}}}{2}+\frac{{\mu}_{\mathrm{IR}}-{\mu}_{\mathrm{IL}}}{2}\mathrm{cos}(\frac{\pi x}{a}\left)\right],$$## (3)

$$\epsilon (\omega ,x)={\epsilon}_{0}{\epsilon}_{R}(\omega )\mathrm{cos}\left(\frac{\pi x}{a}\right)-i{\epsilon}_{0}[\frac{{\epsilon}_{\mathrm{IR}}+{\epsilon}_{\mathrm{IL}}}{2}+\frac{{\epsilon}_{\mathrm{IR}}-{\epsilon}_{\mathrm{IL}}}{2}\mathrm{cos}(\frac{\pi x}{a}\left)\right],$$According to the Kramers-Kronig (KK) relations that are fully valid for metamaterials^{34}^{,}^{35} (i.e., due to causality), the real and the imaginary parts of such complex media are related, but it is an integral interdependence which by no means stipulates that they have to follow the same trends in the relatively narrow range of frequencies typically observed in metamaterials. On the contrary, in a more general case, these parts of the dependencies can assume very different forms, up to the point of assuring a possibility to compensate unavoidable losses in the negative refractive index part through the introduction of active media,^{36} where KK relations still remain valid.

The apparent interdependence between the parameters of the two media is due to the mathematical properties of the transition function model, sinusoidal in this case. However, further away from the transition region, the imaginary parts of the permittivity [${\epsilon}_{\mathrm{IR}}(\omega )$ and ${\epsilon}_{\mathrm{IL}}(\omega )$] and permeability [${\mu}_{\mathrm{IR}}(\omega )$ and ${\mu}_{\mathrm{IL}}(\omega )$] in PRM and NRM media, respectively, can be chosen independently from each other. The reason for giving Eqs. (2) and (3) is just a need to formulate a simple mathematical model of gradual transition between NRM and PRM. In a more elaborate model, which will be the subject of our future work, this distinction between the two materials will become more manifest.

Unlike some other functional dependences of PRM–NRM transitions studied so far, for instance $\mathrm{tan}\text{\hspace{0.17em}}h(x)$ model,^{24} the sinusoidal model provides only for relatively slow transitions and it is not equally obvious where we are “far away from the transition region.” In continuous models, the freedom of choice of parameters e.g., the permittivity ${\epsilon}_{\mathrm{IR}}(\omega )$ and ${\epsilon}_{\mathrm{IL}}(\omega )$ as well as permeability ${\mu}_{\mathrm{IR}}(\omega )$ and ${\mu}_{\mathrm{IL}}(\omega )$ in PRM and NRM media, respectively, is an asymptotic statement. For example, $\mathrm{tan}\text{\hspace{0.17em}}h(\infty )\to 1$ is only an asymptotic constant. The fact is that in a continuous model there is everywhere a spatial dependency and interdependency of material parameters. But at some asymptotic points, we calibrate the spatially constant and frequency-dependent material parameters to correspond to the actual PRM and NRM media far away from the transition region which we desire that they have. For instance, Eqs. (2) and (3) at the maxima of the cosine function (in the middle of PRM), where $\pi x/a=2n\pi $ or $x=2na$ ($n=0,1,2,3,\dots $), give

Thus, we see that in the middle of the PRM, we have the imaginary parts of permittivity and permeability ${\epsilon}_{\mathrm{IR}}(\omega )$ and ${\mu}_{\mathrm{IR}}(\omega )$, respectively, while in the middle of the NRM, we have the imaginary parts of permittivity and permeability ${\epsilon}_{\mathrm{IL}}(\omega )$ and ${\mu}_{\mathrm{IL}}(\omega )$, respectively, and in Eqs. (2) and (3) these parameters can be chosen arbitrarily.

In order to obtain a constant wave impedance throughout the structure, we now introduce a further requirement that the real and imaginary parts of the effective permittivity and permeability satisfy the condition

## (4)

$$\frac{{\mu}_{\mathrm{IR}}(\omega )+{\mu}_{\mathrm{IL}}(\omega )}{2{\mu}_{R}(\omega )-i[{\mu}_{\mathrm{IR}}(\omega )-{\mu}_{\mathrm{IL}}(\omega )]}=\frac{{\epsilon}_{\mathrm{IR}}(\omega )+{\epsilon}_{\mathrm{IL}}(\omega )}{2{\epsilon}_{R}(\omega )-i[{\epsilon}_{\mathrm{IR}}(\omega )-{\epsilon}_{\mathrm{IL}}(\omega )]}=\beta (\omega ),$$The condition (4) is a restrictive mathematical requirement on the complex permittivity and permeability that reduces our analysis to a special case. A justification for the requirement (4) is based on the fact that both permittivity and permeability of many of the NRM structures reported until now can be described by Drude or Lorentz models, i.e., that $\epsilon \left(\omega \right)$ and $\mu \left(\omega \right)$ are strongly resonant and thus quite narrow. In order to obtain the widest possible frequency range of negative refractive index, it is then useful to have the best possible overlap between the ranges of negative values of the real parts of $\epsilon \left(\omega \right)$ and $\mu \left(\omega \right)$. In an ideal situation their dispersions in the resonant range would be therefore identical. On the other hand, in order to preserve causality, the imaginary parts of both $\epsilon \left(\omega \right)$ and $\mu \left(\omega \right)$ must be positive and their dispersions are determined by the real parts—thus the imaginary parts should also overlap. Actually, Eq. (4) can be seen as a condition, although too stringent, for the maximum bandwidth of Drude- or Lorentz-type resonant NRM structures.

Regarding the fabrication of NRM structures with largely overlapping complex $\epsilon \left(\omega \right)$ and $\mu \left(\omega \right)$, it was our reasoning that since NRM are artificial media they allow us, in principle, to design $\epsilon \left(\omega \right)$ and $\mu \left(\omega \right)$ separately, as is the case with prototypical NRM, the split ring resonators combined with the wire media.^{3}^{,}^{9} Using Eq. (4), we readily obtain

## (5)

$$\mu ={\mu}_{0}\frac{{\mu}_{\mathrm{IR}}+{\mu}_{\mathrm{IL}}}{2\beta}\left[\mathrm{cos}\right(\frac{\pi x}{a})-i\beta ],\phantom{\rule[-0.0ex]{2em}{0.0ex}}\epsilon ={\epsilon}_{0}\frac{{\epsilon}_{\mathrm{IR}}+{\epsilon}_{\mathrm{IL}}}{2\beta}\left[\mathrm{cos}\right(\frac{\pi x}{a})-i\beta ].$$Except for condition Eq. (4), our method allows for arbitrary temporal dispersion. Upon condition Eq. (4), the wave impedance $Z={Z}_{0}Z(\omega )=\sqrt{\mu (\omega ,x)/\epsilon (\omega ,x)}$ becomes constant throughout the entire structure; as a result, there is no reflection on the graded interfaces between the two materials. Equation (1) with $\mu =\mu (\omega ,x)$ and $\epsilon =\epsilon (\omega ,x)$ given by Eqs. (2) and (3), respectively, is easily transformed into hypergeometric equations. If such an approach is used, the derivation is rather simple and does not require the use of any specialized software. The solutions to these equations involving the appropriate hypergeometric functions, for the given set of parameters, are reduced to relatively simple elementary functions. For this particular graded index structure, the exact solutions to the two differential equations (2) and (3) reduce to a remarkably simple form

## (6)

$$E(x)={E}_{0}{\mathrm{e}}^{-\kappa \beta x}\mathrm{exp}(-i\frac{\kappa a}{\pi}\mathrm{sin}\frac{\pi x}{a}),\phantom{\rule[-0.0ex]{2em}{0.0ex}}H(x)={H}_{0}{\mathrm{e}}^{-\kappa \beta x}\mathrm{exp}(-i\frac{\kappa a}{\pi}\mathrm{sin}\frac{\pi x}{a}),$$## (7)

$$\kappa =k+i\alpha =\frac{\omega}{c}\sqrt{{\mu}_{R}{\epsilon}_{R}}+i\frac{\omega}{2c}\sqrt{\frac{{\epsilon}_{R}}{{\mu}_{R}}}({\mu}_{\mathrm{IL}}-{\mu}_{\mathrm{IR}}).$$From Eq. (7), one may draw a superficial conclusion that the loss of the entire structure stems from the imaginary part of the permeability of both PRM and NRM only. However, this is not the case, and Eq. (7) is just one way of writing the loss parameters of the structure. Using Eq. (4), it is possible to formulate the equivalent result to Eq. (7) where only the imaginary part of the permittivity of both PRM and NRM appears. Thus, Eq. (7) by no means implies that the contribution from the imaginary part of permittivity is marginal and can be ignored. It is just a matter of an arbitrary choice of presenting the quantities interrelated by Eq. (4).

Although both $\kappa $ and $\beta $ are complex numbers, it should be noted that the product $\kappa \times \beta $ is a real number. From Maxwell’s equation (1), the field amplitudes are related by ${E}_{0}={Z}_{0}Z(\omega ){H}_{0}$. The exact solutions (6) are valid for arbitrary loss factors in NRM and PRM. In the PRM slab around the origin, i.e., in the limit $x\to 0$ we readily obtain the time-domain fields of the form (since $E(x,t)=\mathrm{Re}[E(x){e}^{i\omega t}]$, etc.)

## (8)

$$E(x,t)\sim {E}_{0}{e}^{-{\gamma}_{1}x}\text{\hspace{0.17em}}\mathrm{cos}(\omega t-kx)H(x,t)\sim {H}_{0}{e}^{-{\gamma}_{1}x}\text{\hspace{0.17em}}\mathrm{cos}(\omega t-kx)$$## (9)

$$E(x,t)\sim {E}_{0}{e}^{-{\gamma}_{2}(x-a)}\text{\hspace{0.17em}}\mathrm{cos}[\omega t-(-k)x],\phantom{\rule{0ex}{0ex}}H(x,t)\sim {H}_{0}{e}^{-{\gamma}_{2}(x-a)}\text{\hspace{0.17em}}\mathrm{cos}[\omega t-(-k)x].$$In the results shown in Eqs. (8) and (9) for the fields, the loss factors ${\gamma}_{1}$ and ${\gamma}_{2}$ are given by

## (10)

$${\gamma}_{1}=\frac{\omega}{c}\sqrt{\frac{{\epsilon}_{\mathrm{R}}}{{\mu}_{\mathrm{R}}}}{\mu}_{\mathrm{IR}},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\gamma}_{2}=\frac{\omega}{c}\sqrt{\frac{{\epsilon}_{\mathrm{R}}}{{\mu}_{\mathrm{R}}}}{\mu}_{\mathrm{IL}},$$From the results in Eq. (8), it follows that the asymptotic wavevector in the right-handed material is ${\overrightarrow{k}}_{\mathrm{RHM}}=+k{\overrightarrow{e}}_{x}$, i.e., the wave propagates in the $+x$-direction. On the other hand, from the results (9), it follows that the asymptotic wavevector in the left-handed material is ${\overrightarrow{k}}_{\mathrm{LHM}}=-k{\overrightarrow{e}}_{x}$, i.e., the wave propagates in the $-x$-direction. However, the energy flux (the Poynting vector) is in the $+x$-direction in both media. Thus, in the limit of small $x$, we have the correct wave behavior in both the PRM and the NRM slabs. The wave changes the direction periodically along the periodic structure, which will be obvious from the figures in Sec. 5.

## 4.

## Numerical Model of Metamaterials

A dispersive transmission line matrix (TLM) $Z$-transform model of the lossy MM-composite, described in Ref. 37, is used here to verify the analytical solution for gradient index metamaterials with arbitrary loss factor in PRM and NRM presented in Secs. 2 and 3. This model follows the notation used in Refs. 38 and 39 to describe various types of conventional linear time-dependent materials with the purpose that it can be easily incorporated into the algorithm of the so-called $Z$-transform-based TLM method, given in Ref. 38.

The TLM $Z$-transform model of the lossy MM-composite is based on the Drude dispersive model as it allows to characterize MM-composite response in a much wider frequency range than, e.g., the Lorentz dispersion model. However, it could be easily adapted to describe any higher-order material responses. In this paper, the Drude model describing the frequency dependence of electric and magnetic conductivities is used

## (11)

$${\sigma}_{e}(\omega )=\frac{{\sigma}_{e0}}{1+j\omega {\tau}_{e}}=\frac{{\omega}_{pe}^{2}{\tau}_{e}{\epsilon}_{0}}{1+j\omega {\tau}_{e}},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\sigma}_{m}(\omega )=\frac{{\sigma}_{m0}}{1+j\omega {\tau}_{m}}=\frac{{\omega}_{pm}^{2}{\tau}_{m}{\mu}_{0}}{1+j\omega {\tau}_{m}},$$As an alternative, the Drude model describing the frequency dependence of permittivity and permeability (i.e., electric and magnetic susceptibilities) of MM composites can be used, but both models give identical results as shown in Ref. 37. In addition, using the relations

## (12)

$$\epsilon (\omega )={\epsilon}_{0}[1-j\frac{{\sigma}_{e}(\omega )}{\omega {\epsilon}_{0}}],\phantom{\rule[-0.0ex]{2em}{0.0ex}}\mu (\omega )={\mu}_{0}[1-j\frac{{\sigma}_{m}(\omega )}{\omega {\mu}_{0}}]$$In Ref. 38, the Drude model was also used but only to describe electric conductivity of a nonmagnetized plasma with collisions. Also, instead of the exponential $Z$-transform employed in Ref. 38 to transfer frequency dependence of the considered material properties to the $Z$-domain, the bilinear $Z$-transform is adopted in Ref. 37 to develop TLM $Z$-transform model of lossy MM since the bilinear discretization provides a much more accurate scheme.

In this section, the dispersive TLM $Z$-transform model of the lossy MM composite, based on the Drude model for electric and magnetic conductivities of MM composite, will be briefly described and also presented in a more general form through signal flow diagrams. For the sake of brevity and clarity of the formulation to follow, only diagrams corresponding to the $y$-component of the electric field and the $z$-component of the magnetic field are shown in this paper. The calculation for other field components within the proposed numerical model can also be illustrated using similar diagrams.

In general, the TLM $Z$-transform method algorithm can be described by the following three equations:^{38}

## (13)

$${\underline{F}}^{r}={\underline{\underline{R}}}_{1}^{T}{\underline{V}}^{i}-0.5{\underline{V}}_{f},\phantom{\rule[-0.0ex]{1em}{0.0ex}}\underline{F}=\underline{\underline{t}}(z){\underline{F}}^{r},\phantom{\rule[-0.0ex]{1em}{0.0ex}}{\underline{V}}^{r}=\underline{\underline{R}}\underline{F}-\underline{\underline{P}}{\underline{V}}^{i},$$For the modeling of general linear isotropic frequency-dependent materials, only the transmission block $\underline{\underline{t}}(z)$ needs to be calculated. The second term of Eq. (13) can be then reduced, for the considered field components, to

## (14)

$$4{V}_{y}^{r}=4{V}_{y}+{g}_{e}(z){V}_{y}+4\left(\frac{1-{z}^{-1}}{1+{z}^{-1}}\right){\chi}_{e}(z){V}_{y},$$## (15)

$$-4{i}_{z}^{r}=4{i}_{z}+{r}_{m}(z){i}_{z}+4\left(\frac{1-{z}^{-1}}{1+{z}^{-1}}\right){\chi}_{m}(z){i}_{z}.$$Expressing electric and magnetic conductivities, given by the Drude model in Eqs. (11) and (12), by using the normalized conductivity ${g}_{e}(\omega )={\sigma}_{e}(\omega ){\eta}_{0}\mathrm{\Delta}l$ and normalized resistivity ${r}_{m}(\omega )={\sigma}_{m}(\omega )\mathrm{\Delta}l/{\eta}_{0}$, respectively, and transforming them to the $Z$-domain using the bilinear transformation $j\omega \to 2(1-{z}^{-1})/[\mathrm{\Delta}t(1+{z}^{-1})]$ ($\mathrm{\Delta}l$ and $\mathrm{\Delta}t$ are the space and time discretization steps in TLM, respectively), the following representations in the $Z$-domain are obtained

## (16)

$${g}_{e}(z)=(1+{z}^{-1})\frac{{g}_{ec}}{{B}_{ce}(1-{z}^{-1}{A}_{ce}/{B}_{ce})},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{r}_{m}(z)=(1+{z}^{-1})\frac{{r}_{mc}}{{B}_{cm}(1-{z}^{-1}{A}_{cm}/{B}_{cm})},$$The coefficients shown in Figs. 3 and 4 can be found after arranging that the frequency dependence of electric and magnetic conductivity is represented as a function of the field value at the previous time step by taking partial fraction expansions forms shown below

## (17)

$$(1+{z}^{-1}){g}_{e}(z)={g}_{e0}+{z}^{-1}\lfloor {g}_{e1}+\overline{{g}_{e}(z)}\rfloor ,\phantom{\rule[-0.0ex]{1em}{0.0ex}}(1+{z}^{-1}){r}_{m}(z)={r}_{m0}+{z}^{-1}\lfloor {r}_{m1}+\overline{{r}_{m}(z)}\rfloor .$$Inserting Eq. (16) into Eq. (17) gives

## (18)

$${g}_{e0}=\frac{{g}_{ec}}{{B}_{ce}},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{g}_{e1}=\frac{{g}_{ec}}{{B}_{ce}}(2+{A}_{ce}/{B}_{ce})={g}_{e0}(2+{a}_{ce}),$$## (19)

$$\overline{{g}_{e}(z)}={z}^{-1}\frac{{g}_{ec}(1+2{A}_{ce}/{B}_{ce}+{A}_{ce}^{2}/{B}_{ce}^{2})}{{B}_{ce}(1-{z}^{-1}{A}_{ce}/{B}_{ce})}=\frac{{z}^{-1}{b}_{ce}}{1-{z}^{-1}{a}_{ce}},$$## (20)

$${r}_{m0}=\frac{{r}_{mc}}{{B}_{mc}},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{r}_{m1}=\frac{{r}_{mc}}{{B}_{cm}}(2+{A}_{cm}/{B}_{cm})={r}_{m0}(2+{a}_{cm}),$$## (21)

$$\overline{{r}_{m}(z)}={z}^{-1}\frac{{r}_{mc}(1+2{A}_{cm}/{B}_{cm}+{A}_{cm}^{2}/{B}_{cm}^{2})}{{B}_{cm}(1-{z}^{-1}{A}_{cm}/{B}_{cm})}=\frac{{z}^{-1}{b}_{cm}}{1-{z}^{-1}{a}_{cm}},$$## (22)

$${T}_{ce}={(4+{g}_{e0}+4{\chi}_{e\infty})}^{-1},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{T}_{cm}={(4+{r}_{m0}+4{\chi}_{m\infty})}^{-1},$$## (23)

$${k}_{ce}=-(4+{g}_{e1}-4{\chi}_{e\infty}),\phantom{\rule[-0.0ex]{2em}{0.0ex}}{k}_{cm}=-(4+{r}_{m1}-4{\chi}_{m\infty}).$$The results of the numerical calculations using TLM $\mathrm{Z}$-transform for the present problem are given and compared with the exact analytical results in the next section.

## 5.

## Results and Discussion

The exact analytical solutions for the real part of the electric field $E(x)$, given by Eq. (6), for two different values of the numerical parameters (blue lines) are presented and compared to the corresponding numerical solutions (brown lines) in Figs. 5 and 6.

Both figures show that there is no reflection at the interfaces between the right-handed and the left-handed material slabs in this particular case. This is expected, since in our case, the impedance is constant throughout the entire space. From the presented curves, we see the obvious change of the direction of the wave at the boundaries between the slabs (e.g., at $x=-a/2$ and $x=a/2$). Furthermore, we see in Fig. 5, for ${\gamma}_{2}=7{\gamma}_{1}$ there is considerably stronger attenuation of the signal in NRM compared to that in PRM. On the other hand, in Fig. 6, for ${\gamma}_{2}=1.5{\gamma}_{1}$ the attenuation of the signal in NRM is only slightly stronger compared to that in PRM. Furthermore, both Figs. 5 and 6 show an excellent agreement between the exact analytic results and the corresponding numerical results obtained using TLM $Z$-transform. The obtained numerical and analytical curves are practically indistinguishable.

From the results (6), and utilizing the definition of the Poynting vector,^{40} we readily see that the time average of the Poynting vector through the structure is given by

## (24)

$$\langle \overrightarrow{\mathrm{\Pi}}(x)\rangle =\frac{1}{2}\mathrm{Re}[\overrightarrow{E}(x)\times {\overrightarrow{H}}^{*}(x)]=\frac{1}{2}{E}_{0}{H}_{0}{e}^{-2\kappa \beta x}{\overrightarrow{e}}_{x}\mathrm{.}$$This is an exponentially decaying, strictly positive function, showing that the energy flows in the positive $x$-direction throughout the entire periodical PRM-NRM structure, as expected. In particular, for the case without losses the time average of the Poynting vector is constant throughout the structure and equal to

## 6.

## Conclusion

We have presented a simple exact analytical solution to Helmholtz equation for periodic structures with graded permittivity and permeability profile changing according to a cosine function along the direction of propagation, with arbitrary loss factors in PRM and NRM. We analyzed a special case of matched impedance throughout the structure where the real parts of the effective permittivities and permeabilities have opposite signs and equal absolute values, while the imaginary parts are fully arbitrary. We compared the exact analytical results with the corresponding numerical results obtained using TLM $Z$-transform and obtained an excellent agreement between the analytical and numerical results.

The model is valid for arbitrary temporal dispersion and arbitrary losses as long as the general mathematical and physical constraints are satisfied (e.g., Kramers-Kronig relations).

It should be noted that, throughout the present paper, when we state that the loss factors can be chosen arbitrarily, we tacitly assume that they can be chosen arbitrarily as long as the general mathematical and physical constraints, e.g., Kramers-Kronig relations, are satisfied. Such constraints do not, in general, impose serious restrictions to (theoretically) match the impedance between the two media. The challenge is rather the practical realization of media which satisfy the conditions posed in Eq. (4). The model allows a straightforward generalization to any periodic refractive index profile using suitable Fourier series, including abrupt transitions between two materials, which will be the subject of coming studies. Other future challenges include the study of the extension of the present model to various two-dimensional and three-dimensional structures and the case of arbitrary index profiles.

## Acknowledgments

The work of T.A. and N.D. was funded by the Serbian Ministry of Education, Science and Technological Development within the project TR-32024. The work of Z.J. was funded by the Serbian Ministry of Education, Science and Technological Development within the project TR-32008. M.D. developed the analytical solution and wrote a part of the manuscript. M.N. took part in developing the analytical model. T.A. developed numerical model and performed numerical verification of the results. N.D. developed numerical model and wrote a part of the manuscript. Z.J. took part in developing the concept and wrote a part of the manuscript.

## References

## Biography

**Mariana Dalarsson** is a teaching and research assistant at Royal Institute of Technology, School of Electromagnetic Engineering. She received her MSc in microelectronics and applied physics from the Royal Institute of Technology, Stockholm, Sweden. Her interests include inverse problems in electromagnetics, electromagnetics of stratified media, double-negative metamaterials, mathematical physics, etc. She is an author of 18 peer-reviewed publications, including 10 journal papers and a textbook in theoretical physics.

**Martin Norgren** received his PhD in electrical engineering-electromagnetic theory from the Royal Institute of Technology, Stockholm, Sweden. His interests include electromagnetic theory of stratified media, chiral materials, double-negative metamaterials, etc. He is an author of a number of peer-reviewed publications, including 35 journal papers. He is a professor at Royal Institute of Technology, School of Electromagnetic Engineering.

**Tatjana Asenov** received his Dipl.-Ing. degree from the Faculty of Electronic Engineering, University of Niš, Serbia, in 2008. Since 2009, she has been a research assistant and a PhD candidate at the Department of Telecommunications, Faculty of Electronic Engineering, Serbia. She conducts research, funded by the Serbian Ministry of Education and Science, on metamaterials and computational and applied electromagnetics with strong emphasis on microwave applications.

**Nebojšа Dončov** was born in Niš, Serbia, in 1970. He received his Dipl.-Ing., MSc, and PhD degrees from the Faculty of Electronic Engineering, University of Nis, Serbia, in 1995, 1999, and 2002, respectively. From 1995 to 2001, he was with the Department of Telecommunications, Faculty of Electronic Engineering, Serbia, as a research assistant. From 2001 to 2004, he was working with Flomerics Ltd., Electromagnetics Division, UK, as a research and development engineer. In 2004, he joined the Department of Telecommunications, Faculty of Electronic Engineering, Serbia, where he is now an associate professor. His current research interests are in computational and applied electromagnetics with a particular emphasis on TLM and network methods applications in microwaves and EMC. He was the recipient of the International Union of Radio Science (URSI) Young Scientist Award in 2002.

**Zoran Jakšić** received his Dipl.-Ing., Mag. Sci., and PhD degrees in electrical engineering-engineering physics from the School of Electrical Engineering, University of Belgrade. His interests are nanophotonics and nanoplasmonics, including optical metamaterials, NEMS and MEMS sensors and infrared detectors. He is an author of about 260 peer-reviewed publications including 60 journal papers and five book chapters. He is a full research professor and the science director with the Center of Microelectronic Technologies and Single Crystals, a department of the Institute of Chemistry, Technology and Metallurgy, University of Belgrade, Serbia.