## Abstract

We investigate a novel implementation of hyperbolic metamaterial (HM) at far-infrared frequencies
composed of stacked graphene sheets separated by thin dielectric layers. Using the surface
conductivity model of graphene, we derive the homogenization formula for the multilayer structure by
treating graphene sheets as lumped layers with complex admittances. Homogenization results and
limits are investigated by comparison with a transfer matrix formulation for the HM constituent
layers. We show that infrared iso-frequency wavevector dispersion characteristics of the proposed HM
can be tuned by varying the chemical potential of the graphene sheets via electrostatic biasing.
Accordingly, reflection and transmission properties for a film made of graphene-dielectric
multilayer are tunable at terahertz frequencies, and we investigate the limits in using the
homogenized model compared to the more accurate transfer matrix model. We also propose to use
graphene-based HM as a super absorber for near-fields generated at its surface. The power emitted by
a dipole near the surface of a graphene-based HM is increased dramatically (up to 5 ×
10^{2} at 2 THz), furthermore we show that most of the scattered power is directed into the
HM. The validity and limits of the homogenized HM model are assessed also for near-fields and show
that in certain conditions it overestimates the dipole radiated power into the HM.

© 2013 Optical Society of America

## 1. Introduction

A stack of graphene sheets, separated by subwavelength dielectric spacers, can be regarded as a composite material with uniaxial electric properties under certain conditions. Graphene layers strongly affect the complex effective permittivity of the composite for electric field components polarized parallel to the graphene sheets. Uniaxial anisotropic materials in general offer a variety of interesting electromagnetic properties. In particular, we investigate here a subcategory denoted as hyperbolic metamaterials (HMs), named after the hyperbolic iso-frequency wavevector dispersion that arise due to the negative permittivity experienced by the electric field component along either the axis of anisotropy or a direction orthogonal to the axis of anisotropy [1, 2].

Strong interest in HMs is based on their specific property that enables propagation in a very wide spatial spectrum, that would be otherwise evanescent in free space, which is in principle unbounded for the ideal case of purely hyperbolic iso-frequency wavevector dispersion. In case of realistic hyperbolic-like dispersion, the spatial spectrum allowed for propagation can still be extremely wide, as shown in [3–5]. This property is shared with uniaxial anisotropic materials having elliptic iso-frequency dispersion diagram, with a very large axial ratio. Though in practical cases purely hyperbolic dispersion cannot be obtained, effective medium models based on hyperbolic dispersion proves to be a very useful tool for understanding the physics behind the interesting electromagnetic properties of these metamaterials with extremely subwavelength features.

Recently, metal-dielectric multilayers were proposed as candidates to realize HMs at optical frequencies [3, 6, 7]. The wide spatial spectrum of propagation supported by these HMs can lead to novel phenomena as increasing the power emitted by imposed dipoles [3,8] or scattered by nanoparticles [3,9] on HM surfaces, and this power is mostly directed into HMs. This exotic property enables features like focusing with very subwavelength resolution [10,11], controlling absorption [12], enhancement of spontaneous emission [13], increasing the decay rate of emitters [6], designing quantum and thermal emitters [14]. HMs can also host backwards waves and thus they are utilized for achieving negative refraction as in [15]. As stated in [4], HMs are considered to be promising materials for advancement in the fields of tunable surface plasmon polaritons, super Planckian thermal emission [5], radiative decay engineering [16], and nano-imaging.

HMs attract attention also because they are easy to fabricate using metal-dielectric multilayers or metallic nanowires embedded in dielectric substrates. Also, doped semiconductors [17, 18] and conductive oxides used for generating surface plasmons can be used for HM designs in near- and mid-IR frequency bands [15]. In multilayer structures, metal is used as a negative permittivity layer-spaced by dielectric layers, overall creating a negative permittivity effect for the electric field components tangential to the layers. This effect does not rely on any resonant behavior and thus is a very wide-frequency-band.

In this paper we investigate a multilayer HM design based on stacking dielectric layers and graphene sheets. A graphene sheet (a one-atom-thick carbon layer) is able to support surface plasmon modes at terahertz frequencies [19–23]. A recent work in [24] a discussed HM based on a semi-infinite stack of graphene-dielectric multilayers, studied at a temperature of 4 K, thus assuming graphene losses negligible, leading to purely real permittivity and wavenumbers. Our analysis instead accounts for losses at a room temperature of 300 K and their effect on (i) the effective permittivity that assumes complex values, (ii) the hyperbolic-like dispersion (it is not exactly hyperbolic), hence (iii) wavenumbers that may assume complex values, and (iv) on the radiated power by a dipole near the HM surface, where losses play an important role. Moreover, we investigate practical cases with a finite number of graphene-dielectric layers, and quantitatively show tunability aspects of graphene-based HMs using electrostatic biasing. We also show a detailed study about the dependence of power spectrum emitted by a dipole source in the proximity of the graphene-based HM on the number of layers, as well as on frequency. Alternative to the HM implemented as a graphene-dielectric multilayer studied here, in [25] graphene stripes analogous to a metallic wire medium, are utilized for realizing hyperbolic dispersion in cylindrical wavenumber coordinates, with the aim of designing a hyperlens. The peculiar electronic properties of graphene [26, 27], have been investigated for different electromagnetic applications such as lensing [28], transformation optics [29], nanomechanical resonators [30], and solar cells [31]. Moreover, based upon periodic patterning of graphene, bi-periodic and/or multilayered graphene structures were extensively studied for enhanced transmission [32], optical absorption [33], and tunable metasurfaces in [34], as well as isolators and polarizers in the microwave regime [35–37]. It was shown that crystalline Graphite (the 3D parent of graphene) possess indefinite permittivity at UV frequencies [38].

A graphene sheet has properties at THz frequencies similar to those of a thin metal film at optical frequencies, as shown in our work. In principle, any inductive infinitesimally-thin layer can be used to realize HMs, however, in the terahertz regime, designing inductive layers is difficult due to metallic losses and spatial and frequency dispersion introduced by periodically patterned conductive layers. For this aim, metallic meshes are mainly effective only in lower microwave regime [39–41]. On the other hand, the use of highly dispersive metals is practical the optical frequencies below the plasma frequency. Here we show that stacking graphene sheets can be utilized for designing tunable HMs in a wide frequency spectrum ranging from millimeter-waves up to tens of terahertz frequencies, encompassing the whole far-infrared band. In Section 2, we develop an effective medium approximation (EMA) to facilitate the characterization of graphene-based multilayer structure and its use in possible devices and we assess its limits. Moreover, in Section 2 we describe the basic properties of graphene and the associated hyperbolic wavenumber dispersion of the multilayer structure composed of graphene-dielectric layers. In Section 3 we investigate plane wave transmission and reflection for a thin slab made of several graphene sheets, as well as their tunability features, and we also show how EMA is able to describe these properties. In Section 4 we investigate the radiation of a dipole at the interface of a finite thickness HM and we show that the HM is able to enhance the total power radiated by several orders of magnitude, reporting an enhancement of 5 × 10^{2} at 2 THz. We also show that most of the power is directed into the HM, offering a viable route for wide band and wide incidence-angle super absorption interfaces at far-infrared frequencies, as previously discussed in [3, 9, 13, 42] for optical frequencies. We show that this large enhancement of power emission is associated to the wide spatial spectrum being able to propagate inside the HM, that would be otherwise evanescent in free space.

## 2. Tunable hyperbolic metamaterial made of graphene-dielectric multilayers

A graphene monolayer is electrically characterized by its surface conductivity *σ*(*ω*, *μ _{c}*), where

*μ*is the chemical potential related to the electrostatic biasing, which quantifies the electronic transport properties [43]. The frequency dependent conductivity follows the interband (bound-electrons) and the intraband (free-electrons) sum rules [44, 45]. Spatial dispersion of graphene has negligible effects, since the graphene lattice constant

_{c}*a*≈0.246 nm is extremely subwavelength at THz frequencies [34]. Therefore, a graphene layer is modeled by the local isotropic surface conductivity

*σ*(

*ω*,

*μ*) =

_{c}*σ*′ +

*jσ*″ (assuming time-harmonic fields of the form

*e*) that is calculated by the Kubo formula

^{jωt}*k*is the Boltzmann constant,

_{B}*h*is the Planck constant,

*f*(

_{D}*ζ*) = [

*e*

^{(ζ−μc)/(kBT)}+ 1]

^{−1}is the Fermi-Dirac distribution, and Γ is the phenomenological scattering rate. Throughout our discussion, we assume Γ = 0.1 meV, which is within the range of values considered in other studies [21, 37, 46, 47], at room temperature

*T*= 300 K. Consider a periodic stack of graphene-dielectric layers, as depicted in Fig. 1. Each dielectric spacer has a subwavelength thickness

*d*and relative permittivity

*ε*. In this premise, we assume that graphene sheets are electronically isolated, i.e., the electronic band structure of a graphene sheet is not affected by the neighboring graphene sheets (interlayer electronic coupling mechanisms as well as tunneling effects are ignored, due to the significant thickness of the dielectric with respect to quantum scales). We model graphene sheets as lumped complex-admittance layers, due to their extremely sub-wavelength thickness. Wave propagation in the multilayer structure, depicted in Fig. 1, can be modeled using (i) EMA that models the multilayer as a homogeneous medium and (ii) by applying Bloch theory (as in Chapter 8 in [48]) using the transfer matrix of a unit cell. When applying EMA, the multilayered structure is modeled as a homogeneous uniaxial anisotropic medium (with axis of anisotropy along

_{d}*z*) whose effective relative permittivity tensor

*ε̱*_{eff}is a diagonal matrix in Cartesian coordinates given as

Since a graphene sheet is infinitesimally thin with respect to the dielectric thickness, one may assume that *ε _{z}* =

*ε*, because

_{d}*z*–directed electric field would not excite any current in the graphene sheet. The transverse effective relative permittivity

*ε*is determined as follows, considering a unit cell made of a dielectric layer with dielectric constant

_{t}*ε*between

_{d}*z*= 0 and

*z*=

*d*, and a graphene sheet at

*z*= 0. Within this unit cell one can write ∇ ×

**H**=

*jωε*

_{0}

*ε*

_{d}**E**+

**J**=

*jωε*

_{0}

*ε̱**·*

_{d}**E**, where the current density

**J**[A/m

^{2}] in the graphene sheet is reduced to the surface current along the sheet

**J**=

*δ*(

*z*)

*σ*

**E**

*, and*

_{t}**E**

*is the transverse component of the electric field. Therefore, one has ${\underset{\_}{\mathit{\epsilon}}}_{d}={\epsilon}_{d}\underset{\_}{\mathbf{I}}-j\frac{\sigma}{\omega {\epsilon}_{0}}\delta \left(z\right)\left(\widehat{\mathbf{x}}\widehat{\mathbf{x}}+\widehat{\mathbf{y}}\widehat{\mathbf{y}}\right)$, where*

_{t}**I̱**is the unit tensor, that when averaged over a period along

*z*leads to the effective relative “transverse” permittivity

*ε*,

_{t}The formula for *ε _{t}* in Eq. (3) could be obtained alternatively by following the method used for homogenization involving thin metal-dielectric layers [49]. Accordingly, a graphene sheet may be treated as a layer with extremely subwavelength, but finite, thickness with bulk properties. Exploiting the continuity of the electric field, along

*x*and

*y*, at the boundaries between graphene and dielectric layers, and averaging the transverse component of the effective displacement current over a period

*d*also leads to Eq. (3) (here the effective displacement current account for the displacement current in the dielectric and the conduction current in the graphene sheets). In Eq. (3), we have highlighted that the graphene conductivity is strongly dependent on the frequency and chemical potential. It is important to note that if a graphene sheet has a sufficiently large inductive susceptance, i.e., if

*σ*″ < −

*ωε*

_{0}

*ε*, then the effective relative permittivity term

_{d}d*ε*=

_{t}*ε*′

*−*

_{t}*jε*″

*has a negative real part, i.e.,*

_{t}*ε*′

*< 0. Under this condition, and recalling that*

_{t}*ε*> 0, extraordinary waves, with TM

_{z}*polarization (magnetic field transverse to*

^{z}*z*), are allowed to propagate inside the HM, with wavevectors exhibiting iso-frequency hyperbolic dispersion as explained in [3]; whereas ordinary waves with TE

*polarization (electric field transverse to*

^{z}*z*) are mainly evanescent. This allows for the propagation of TM

*waves with transverse wavenumber ${k}_{t}>\omega \sqrt{{\mu}_{0}{\epsilon}_{0}{\epsilon}_{d}}$, that would be otherwise evanescent in free space.*

^{z}To better illustrate the possible homogenized parameters that can be obtained, we report in Fig. 2 the real and imaginary parts of *ε _{t}* =

*ε*′

*−*

_{t}*jε*″

*using the EMA formula Eq. (3), assuming that graphene sheets are spaced by silica layers with permittivity*

_{t}*ε*= 2.2 and thickness

_{d}*d*= 0.1 μm. Note that

*ε*′

*approaches*

_{t}*ε*at higher frequencies due to saturation of the graphene conductivity to its universal value

_{d}*πe*

^{2}/(2

*h*) [44] which gives a negligible contribution compared to

*ε*. Moreover, the imaginary part of the effective permittivity term (

_{d}*ε*″

*) is remarkably small in a certain frequency band, showing that waves in this composite material can propagate with minimal attenuation. Indeed, it is known that a graphene layer may support weakly attenuated plasmonic modes in at terahertz frequencies [45]. We note that*

_{t}*ε*is very sensitive to the chemical potential

_{t}*μ*, and we show that the zero-crossing frequency of

_{c}*ε*′

*, occurring when*

_{t}*σ*″ = −

*ωε*

_{0}

*ε*, occurs at 6.6, 11, and 24.6 THz, for different values of

_{d}d*μ*= 0, 0.1, and 0.4 eV, respectively. The tunability of the proposed HM structure is then quantified, as seen in Figs. 3(a) and 3(b), where the chemical potential and stacking density (thickness of dielectric layer

_{c}*d*) are varied assuming a frequency of 12 THz. When the dielectric thickness

*d*is increased, one should note that

*ε*′

*increases toward 0. Note also that as the chemical potential (tuned by electrostatic biasing) increases*

_{t}*ε*′

*takes smaller values. Hence, at 12 THz,*

_{t}*ε*′

*is positive (≈2.2) for zero chemical potential, and at higher chemical potential values,*

_{t}*ε*′

*becomes negative: for example at*

_{t}*μ*=0.5 eV, with

_{c}*d*= 0.1 μm, one has

*ε*′

*≈−9.5. In summary, a composite material made by layered graphene sheets possesses desirable performance in terms of losses, inductive response, and tunability from millimeter-waves up to mid-infrared frequencies, hard to find in any other known material. This makes it a good candidate for realizing HM designs in the THz range.*

_{t}When considering plane waves propagating in a homogeneous uniaxial anisotropic medium it is useful to decompose them into the modal polarizations TE* ^{z}* and TM

*. A description using EMA, when valid, gives a neat physical insight into wave propagation in this structure. Wavevector dispersion relations for ordinary (TE*

^{z}*) and extraordinary (TM*

^{z}*) waves in a uniaxial anisotropic materials are written as*

^{z}*z*axis. It is apparent that TM

*waves in a medium with*

^{z}*ε*< 0 exhibit an iso-frequency wavenumber dispersion with hyperbolic shape, as explained in [3,5,50]. This allows the propagation of the extraordinary waves (TM

_{t}*) with any transverse wavenumber ${k}_{t}>\sqrt{{\epsilon}_{d}}{k}_{0}$, that would be otherwise evanescent in a homogeneous dielectric with permittivity*

^{z}*ε*. However, one should note that TE

_{d}*waves are mainly evanescent when*

^{z}*ε*< 0, for any

_{t}*k*. The proposed graphene-dielectric metamaterial can be used for realizing hyperbolic dispersion with

_{t}*ε*< 0 and

_{t}*ε*> 0, and with the present implementation it is not possible to have anisotropy such that

_{z}*ε*> 0 and

_{t}*ε*< 0. Elliptic dispersion regime occurring when both

_{z}*ε*> 0 and

_{t}*ε*> 0, inherently implies a propagating spectrum with ${k}_{t}<\sqrt{{\epsilon}_{z}}{k}_{0}$. Considering the multilayered structure depicted in Fig. 1, we shall consider the root of ${k}_{z}^{2}$, solution of Eq. (5), that corresponds to a wave whose Poynting vector is directed toward the graphene-based HM, i.e, in the −

_{z}*z*direction, as shown with

**v**

*in Fig. 1. We assume here that the*

_{g}*z*-directed wavenumber may assume complex values, i.e.,

*k*=

_{z}*β*−

_{z}*jα*, since the graphene conductivity

_{z}*σ*is complex, modeling the inhomogeneous plane wave spectrum. Accordingly, a wave that carries power in the −

*z*direction shall have the attenuation constant (

*α*) with negative sign, associated to field decay (due to possible losses) along the −

_{z}*z*direction. In general

*β*can have both signs, though in our case it is positive, implying that the TM

_{z}*mode is a backward wave for ${k}_{t}<\sqrt{{\epsilon}_{d}}{k}_{0}$, since*

^{z}*β*< 0 [51, 52].

_{z}α_{z}With the aim of assessing the validity of EMA in Eq. (4), we calculate the iso-frequency wavevector dispersion with the more accurate Bloch theory for the periodic structures. This is done by treating each graphene layer as a complex lumped admittance *Y _{s}* =

*σ*(where the subscript “s” denotes surface) as a shunt load in a transverse equivalent network (TEN, see Chapters 2 and 3 in [1]). This leads to the dispersion relation as (see Appendix A for more details)

*z*-directed wavenumber of a wave inside the dielectric, ${Z}_{d}^{\text{TE}}=\omega {\mu}_{0}/{\kappa}_{d}$, and ${Z}_{d}^{\text{TM}}={\kappa}_{d}/\left({\omega}_{0}{\epsilon}_{0}{\epsilon}_{d}\right)$ are the characteristic wave impedances for TE

*and TM*

^{z}*waves, respectively. Here we report only the dispersion curves that belong to TM*

^{z}*modes for brevity, since they are those exhibiting hyperbolic-like dispersion.. We report in Figs. 4(a)–4(b), and Figs. 4(c)–4(d) the plots of*

^{z}*k*=

_{z}*β*−

_{z}*jα*versus

_{z}*k*, for TM

_{t}*waves at 2 and 12 THz, respectively, by applying EMA Eq. (5) and Bloch theory Eq. (6), for the graphene-dielectic HM (with*

^{z}*ε*= 2.2 and

_{d}*d*= 0.1 μm), for various chemical potentials

*μ*. We plot only the dispersion branch relative to power propagation in the downward direction (see Fig. 1). However one should note the ±

_{c}*k*symmetry in the solutions of Eq. (6). One can observe in Fig. 4 that EMA and Bloch theory are in very good agreement for a wide range of transverse wavenumber

_{z}*k*showing a hyperbolic relation, whereas the curves obtained from the two methods diverge for large

_{t}*k*and the dispersion curve obtained via Bloch theory shows a switching to a mainly evanescent spectrum after certain

_{t}*k*. This observation is in accordance with the simplification of the dispersion relation obtained from Bloch theory as follows. Let us consider the special but important case with |

_{t}*κ*| ≪ 1, and |

_{d}d*k*| ≪ 1, i.e., the period

_{z}d*d*is subwavelength, with respect to the wavenumber in the dielectric and with respect to the Bloch wavenumber (the second inequality also implies that

*k*is far from the edge of the first Brillouin zone where

_{z}*β*= ±

_{z}*π*/

*d*). Thus, we approximate Eq. (6) by taking into account the first and second order Taylor expansion terms corresponding to the approximations cos(

*κ*) ≈ 1 − (

_{d}d*κ*)

_{d}d^{2}/2 and sin(

*κ*) ≈

_{d}d*κ*, that lead to

_{d}dAfter substituting the characteristic impedance by its corresponding value for TE* ^{z}* and TM

*waves, Eq. (7) leads to*

^{z}

^{z}

^{z}It is clear from the analytical analysis above and from Fig. 4 that EMA well describes the wavevector dispersion when |*κ _{d}d*| ≪ 1 and |

*k*| ≪ 1, which can be verified when the transverse wavenumber

_{z}d*k*is not too large. For larger and larger values of

_{t}*k*, the two assumptions would not be valid anymore. In Fig. 4, we also report the evolution of the dispersion curves by varying the chemical potential

_{t}*μ*at 2 THz and 12 THz. For example at 2 THz, when

_{c}*μ*= 0.4 eV, the spatial spectrum is bounded by

_{c}*k*≈ 36

_{t}*k*

_{0}after which

*α*increases dramatically. However, by increasing the chemical potential,

_{z}*ε*′

*assumes larger negative values and*

_{t}*β*−

_{z}*k*dispersion in Fig. 4 evolves to a flatter curve, thus the Brillouin zone edge is reached at larger values of

_{t}*k*. Note that at larger spatial spectrum, the attenuation constant

_{t}*α*increases due to finite losses, as shown in Fig. 4(b). By tuning the chemical potential, the dispersion characteristics can be controlled, for example, at 12 THz the dispersion for unbiased graphene is elliptic as well as when

_{z}*μ*= 0.1 eV. This behavior appears since

_{c}*ε*′

*exhibits zero crossing and becomes positive at 6.6 THz and 11 THz, when*

_{t}*μ*= 0 and

_{c}*μ*= 0.1 eV, respectively. However, when the chemical potential is increased to 0.4 eV, hyperbolic dispersion arises at 12 THz, as shown in Figs. 4(c) and 4(d) where the inset of Fig. 4(c) shows the elliptic behavior for ${k}_{t}<\sqrt{{\epsilon}_{d}}{k}_{0}$. At high frequencies where the conductivity saturates to its universal value, the TM

_{c}*plasmonic modes are extremely confined to graphene layers (*

^{z}*σ*″ becomes very small) and higher frequencies, once

*σ*″ > 0, graphene layers are incapable of supporting those modes [21]. Hence, the spectrum ${k}_{t}>\sqrt{{\epsilon}_{d}}{k}_{0}$ becomes mainly evanescent at frequencies with

*ε*′

*> 0. In other words, after*

_{t}*ε*′

*exhibits a zero-crossing and becomes positive, the wavevector dispersion becomes elliptic.*

_{t}## 3. Plane wave reflection and transmission by a thin film made of graphene-dielectric multilayers

A finite thickness graphene-dielectric multilayer film is considered comprising *N* graphene sheets stacked with silica SiO_{2} dielectric spacer, such that a graphene sheet is at the topmost layer. The thickness of each SiO_{2} spacer is 0.1 μm, and the total multilayered film thickness is *D* = *Nd*. For simplicity, all graphene layers are biased equally using a constant electrostatic potential [34]. For practical consideration, suppose that a thin film of silica (in the order of 100 nm) is deposited on an epitaxially-grown graphene monolayer repeatedly until creating an *N* layer stack; though larger thicknesses could be considered, it is rather simple to achieve the biasing range (*μ _{c}* up to 0.5 eV) using relatively lower electrostatic potential for smaller thicknesses [34, 53].

We investigate reflection and transmission under normal plane wave incidence, and at 30° oblique incidence for both TE* ^{z}* and TM

*polarizations (here*

^{z}*k*=

_{t}*k*

_{0}sin(

*θ*), where

^{i}*θ*is the incidence angle). Reflection and transmission coefficients are reported using the transfer matrix method (solid lines), and using EMA (circles) as well, for the two cases with

^{i}*N*= 10 (

*D*= 1 μm) and

*N*= 20 (

*D*= 2 μm), when

*μ*= 0 eV and 0.4 eV in Fig. 5 and Fig. 6, respectively. As discussed in Section 2, at lower frequencies, TE

_{c}*waves are evanescent for any*

^{z}*k*when

_{t}*ε*′

*< 0, while for TM*

_{t}*waves the iso-frequency wavenumber dispersion is hyperbolic (when*

^{z}*ε*′

*< 0), consequently, a plane wave impinging on the structure with ${k}_{t}<\sqrt{{\epsilon}_{d}}{k}_{0}$ is very weakly transferred, specifically by evanescent coupling. This property is demonstrated in Fig. 5 at frequencies lower than 6 THz, and in Fig. 6 at frequencies lower than 24.6 THz. However, after*

_{t}*ε*′

*exhibits a zero-crossing and becomes positive, the plane wave is able to propagate, and the hyperbolic wavevector dispersion turns into elliptic (thus waves with ${k}_{t}<\sqrt{{\epsilon}_{d}}{k}_{0}$ can propagate). The transmission peak for normal incidence occurs when the effective*

_{t}*ε*′

*is near unity (matched to the free space, where the losses are negligible) for*

_{t}*μ*= 0 eV at ≈ 8 THz (reported in Fig. 5) and for

_{c}*μ*= 0.4 eV at ≈ 33 THz (reported in Fig. 6); this is in accordance with the effective

_{c}*ε*′

*plotted in Fig. 2(a). It is clear that changing the chemical potential of the graphene layers offers great tunability and possibility to control the transmission peak and spectrum. EMA is a good tool to describe plane wave interaction with a graphene-dielectric multilayer thin film, for small dielectric thickness*

_{t}*d*. In order to explore the validity of EMA for thicker dielectric spacers, we report in Fig. 7 the reflection and transmission coefficients for 10 layers of graphene-dielectric layers with varying thickness

*d*, at 10 THz, assuming

*μ*= 0.4 eV. It is shown that EMA yields a noticeable deviation from transfer matrix calculations when

_{c}*d*> 0.2

*λ*

_{0}. Note that in Fig. 7 the transition from hyperbolic to elliptic dispersion occurs at

*d*= 0.02

*λ*

_{0}, which implies

*ε*′

*≈ 0.*

_{t}In particular, these results show two remarkable facts: (i) EMA agrees well with transfer matrix calculations for a wide range of frequencies and dielectric thicknesses, (ii) transmission and reflection by the graphene-based multilayered structure can be effectively tuned by electrostatic biasing. It is evident that graphene layers despite controlling transmission with such small thicknesses, at the same time can be designed to be almost transparent to plane wave excitation [32].

In this Section we have investigated reflection and transmission, and how this is predicted by EMA, for an incident plane wave, however a source or scatterer near the HM interface is able to generate a very wide spatial spectrum of plane waves, including the spectrum with *k _{t}* >

*k*

_{0}, which would be evanescent in free space. In the next Section we show how this wide spectrum is able to propagate inside the HM, similarly to what was done in [3] for a HM at optical frequencies made of dielectric and metallic layers.

## 4. Enhancement of emitted power by an impressed dipole at the surface of a graphene-based HM film

We investigate the power emitted by a transverse dipole located at a distance *h* above the graphene-silica multilayered HM as depicted in Fig. 8(a), over a a silicon substrate (sufficiently thick to be assumed infinitely long, with relative permittivity *ε*_{Si}=11.7). We assume here a unit cell of the HM consisting of a 0.1-μm thick silica layer stacked with a sheet of graphene sheet on top. We calculate the power emitted by the transverse dipole located at *z* = 0 as in Fig. 8(a), by using the spatial spectral formalism of TE* ^{z}* and TM

*waves as outlined in [1]. The total power*

^{z}*P*

_{tot}=

*P*

_{up}+

*P*

_{down}emitted by the transverse dipole illustrated in Fig. 8 is decomposed into the power terms directed toward the +

*z*and −

*z*directions (

*P*

_{up}and

*P*

_{down}, respectively) that are found by the spectral integrals

*p*(

*k*) is the spectral power either in the “up” or “down” direction, with transverse wavenumber

_{t}*k*, and |

_{t}**p**

*| is magnitude of the transverse electric dipole moment. The spectral power ${p}_{\text{up},\text{down}}^{\text{TE},\text{TM}}\left({k}_{t}\right)$ can be written as*

_{t}*Y*represents the equivalent admittance of TE

*/TM*

^{z}*waves seen at the location of the dipole either toward free space or toward the HM (indicated by the subscripts “up” and “down”, respectively), whereas*

^{z}*Y*

_{tot}=

*Y*

_{down}+

*Y*

_{up}, and “*” denotes the complex conjugate. In particular, for what concerns free space (up), the terms ${Y}_{\text{up}}^{\text{TE}/\text{TM}}$ are straightforwardly the TE

*and TM*

^{z}*wave admittances in free space give by ${Y}_{\text{up}}^{\text{TE}}={Y}_{0}^{\text{TE}}={\kappa}_{0}/\left(\omega {\mu}_{0}\right)$, and ${Y}_{\text{up}}^{\text{TM}}={Y}_{0}^{\text{TM}}=\left(\omega {\epsilon}_{0}\right)/{\kappa}_{0}$ where ${\kappa}_{0}=\sqrt{{k}_{0}^{2}-{k}_{t}^{2}}$ is the wavenumber along the*

^{z}*z*axis in free space. The calculation of the admittance

*Y*

_{down}, for either TE

*or TM*

^{z}*waves, is done by translating*

^{z}*Y*

_{HM,}

*, which is the admittance toward −*

_{N}*z*direction, shown in Fig. 8(b), evaluated at the surface of HM (at

*z*= −

*h*), to

*z*= 0 by the simple formula

*Y*

_{HM,N}(see Fig. 8), for either TE

*or TM*

^{z}*waves, is done by using the transfer matrix of*

^{z}*N*unit cells, and representing the silicon substrate at the bottom with a TE

*/TM*

^{z}*wave admittance, as detailed in Appendix B. When using EMA, the multilayer structure is treated as an anisotropic dielectric with relative permittivity in Eq. (2). In Fig. 9 and Fig. 10 (for*

^{z}*μ*= 0 eV and

_{c}*μ*= 0.4 eV, respectively, and assuming a dipole distance

_{c}*h*= 2 μm) we report two power ratios aiming at showing their enhancement: (i) the total power

*P*

_{tot}=

*P*

_{up}+

*P*

_{down}emitted by the dipole normalized by the power emitted by the same dipole in free space

*P*

_{free space}; (ii) the ratio of the power directed downward to the HM,

*P*

_{down}, and the power directed into the upper free space,

*P*

_{up}, for four different cases where the number of graphene sheets is changed as

*N*= 1, 3, 10 and

*N*→ ∞, as well as for a transverse dipole at a distance

*h*above a silicon substrate (dashed lines) for comparison purposes.

The ratio *P*_{tot}/*P*_{free space} also represents the increase of the local density of states (LDOS) with respect to LDOS at a point in free space [5, 50, 54], and this is also referred to as Purcell effect [6, 24]. We model the HM with thickness *Nd* via both the more accurate transfer matrix method (denoted by lines in the figure) and EMA (denoted by circles), and provide the results in Fig. 9(a) and (b) for unbiased graphene (*μ _{c}*= 0 eV), and in Fig. 10 for biased graphene (

*μ*= 0.4 eV). In Fig. 9 one can observe that, at the lowest frequency 0.1 THz, there is a clear trend showing that when the number of layers (

_{c}*N*) increases the ratio

*P*

_{down}/

*P*

_{up}also increases from ≈ 10

^{6}up to ≈ 2.7 × 10

^{8}, when using calculations based on the transfer matrix method. Moreover in the lower frequency range, EMA overestimates

*P*

_{down}/

*P*

_{up}by almost one order of magnitude, as also discussed in [3] for a different HM configuration; however, as the frequency increases EMA and the transfer matrix method agree well. In Fig. 9(b), we observe the same disagreement of the transfer matrix calculations and EMA at lower frequencies, and it is clearly seen that the normalized emitted power (

*P*

_{tot}/

*P*

_{free space}) is ≈ 4 × 10

^{3}at the lowest frequency 0.1 THz and drops linearly as the frequency increases.

*P*

_{down}/

*P*

_{up}exhibits a very sharp drop for

*N*= 1 after around 1 THz, for

*N*= 3 after about 2 THz, for

*N*= 10 after 4 THz, and for

*N*→ ∞ after 6 THz (note that the hyperbolic to elliptic dispersion curve transformation occurs at 6.6 THz when

*μ*= 0 eV obtained via EMA, see Fig. 2, in very good agreement with the

_{c}*N*→ ∞ case). When the chemical potential is increased to

*μ*= 0.4 eV, the “transverse” permittivity

_{c}*ε*′

*decreases to further negative values and the frequency of hyperbolic to elliptic dispersion curve transformation shifts from 6.6 THz (*

_{t}*μ*= 0 eV) to 24.6 THz (

_{c}*μ*= 0.4 eV). In Fig. 10(a), at lower frequencies one can observe that

_{c}*P*

_{down}/

*P*

_{up}is increased by one order of magnitude for all cases whereas

*P*

_{tot}/

*P*

_{free space}decreases by almost one order of magnitude when compared to the case with

*μ*= 0 eV. Moreover the frequency where

_{c}*P*

_{tot}/

*P*

_{free space}exhibits a sharp decrease shifts to a higher frequency when the chemical potential is increased to 0.4 eV in agreement with the change in the frequency of hyperbolic to elliptic dispersion curve transformation, as illustrated from Fig. 4(c). As seen in Figs. 10(a) and 10(b) both the enhancement of emitted power and the ratio of power directed to the –z direction are much larger than the Si-substrate case for a wide frequency band (1–6 THz) in the presence of graphene-dielectric. The interesting features in Figs. 9–10 are related to the power spectrum in Eq. (15), which is described in the following.

We report the emitted power spectrum for TM* ^{z}* waves (solid lines),
${p}^{\text{TM}}\left({k}_{t}\right)={p}_{\text{up}}^{\text{TM}}\left({k}_{t}\right)+{p}_{\text{down}}^{\text{TM}}\left({k}_{t}\right)$, versus normalized transverse wavenumber

*k*/

_{t}d*π*at 0.1 in Fig. 11(a) and 3 THz in Fig. 11(b), varying the number of graphene-dielectric layers assuming

*μ*= 0 eV, and for a better visualization we provide Figs. 11(a) and 11(b) in both logarithmic and linear scales for the horizontal axis, in the left and right panels, respectively. For comparison we also show the power spectrum ${p}^{\text{TE}}\left({k}_{t}\right)={p}_{\text{up}}^{\text{TE}}\left({k}_{t}\right)+{p}_{\text{down}}^{\text{TE}}\left({k}_{t}\right)$ for

_{c}*N*=1 and

*N*→ ∞ (dashed lines). At these two frequencies the composite multilayer exhibits hyperbolic dispersion for TM

*waves and propagation inside the HM occurs for ${k}_{t}>\sqrt{{\epsilon}_{d}}{k}_{0}$. We observe in Fig. 11(b) that in the high*

^{z}*k*spectrum, there are a larger number of spectral peaks when

_{t}*N*increases; eventually yielding a continuous distribution of large spectral intensities when

*N*→ ∞. This explain the advantage of having a large number of layers. Moreover when

*N*→ ∞, one can notice that the power spectrum starts to rise strongly after ${k}_{t}=\sqrt{{\epsilon}_{d}}{k}_{0}$, in agreement with the propagating spectrum’s lower limit in the hyperbolic dispersion diagrams in Fig. 4. We would like to emphasize that plots in linear scale in Fig. 11 clearly show that the propagating power spectrum in the large

*k*region is very wide and therefore it strongly contributes to the spectral integral in Eq. (14). Note that the upper limit of the power spectrum, cut-off at

_{t}*k*

_{t,max}, can be determined by evaluating the evanescent decay in free space between the dipole at

*z*= 0 and the surface of the composite material at

*z*= −

*h*, given by $\text{exp}\left(-\sqrt{{k}_{t}^{2}-{k}_{0}^{2}}h\right)$. For example, by setting $\text{exp}\left(-\sqrt{{k}_{t,\text{max}}^{2}-{k}_{0}^{2}}h\right)=\xi $, where

*ξ*is a predetermined small number, we can consider the power spectrum negligible when

*k*>

_{t}*k*

_{t,max}. It is important to note that for

*h*≪

*λ*

_{0}this upper wavenumber limit

*k*

_{t,max}is independent of the operating frequency when

*k*

_{t,max}≫

*k*

_{0}, because in this case $\text{exp}\left(-\sqrt{{k}_{t,\text{max}}^{2}-{k}_{0}^{2}}h\right)\approx \text{exp}\left(-{k}_{t,\text{max}}h\right)$. These considerations explain why all spectral curves decay with very similar profile for very large

*k*and therefore it is imposed mainly by the dipole distance

_{t}*h*, for both frequencies examined in Fig. 11(a) and (b), i.e., at 0.1 and 3 THz.

At low frequency, in Fig. 11(a) all curves with different *N*, tend to exhibit the same behavior at the large *k _{t}*, in particular when (

*k*/

_{t}d*π*) > 10

^{−2}. The reason of this low frequency property, that does not occur at higher frequency in Fig. 11(b), is explained as follows. Propagation inside such multilayer stack consists in strong evanescent coupling between adjacent graphene sheets [23, 32], that is approximately proportional to $\text{exp}\left(-\sqrt{{k}_{t}^{2}-{\epsilon}_{d}{k}_{0}^{2}}d\right)$. Considering now the large

*k*region where the power spectrum

_{t}*p*

^{TM}(

*k*) is intense in Fig. 11(a), the exponential interlayer decay becomes stronger, at fixed

_{t}*k*, when the frequency decreases (i.e., when ${\epsilon}_{d}{k}_{0}^{2}$ decreases). Thus, at low frequencies, less power is coupled to lower graphene layers, whereas most of the power is coupled to losses in the first graphene sheet closest to the dipole, implying that the number of layers becomes less effective on the power spectrum for large

_{t}*k*and hence on the total emitted power integral in Eq. (14). For example, the case with 0.1 THz in Fig. 11(a), the total TM

_{t}*power emitted by the dipole is dominated by the wide power spectrum region (5 × 10*

^{z}^{−3}

*π*/

*d*) <

*k*<

_{t}*k*

_{t,max}, which is weakly dependent on the number of layers. Under this low frequency condition, we observe that the total TM

*emitted power becomes proportional to*

^{z}*ω*

^{2}, independently on the number of layers, in agreement with the findings in [21] for a single graphene layer. A similar trend occurs for the power emitted as TE

*waves at these lower frequencies, though it is several orders of magnitude weaker than TM*

^{z}*cases for large*

^{z}*k*(see the dashed curves in Fig. 11). Note that, instead, the free space emitted power by a dipole |

_{t}**p**

*|*

_{t}^{2}

*ω*

^{4}

*η*

_{0}/ (12

*πc*

^{2}) is proportional to

*ω*

^{4}. Comparing the low frequency trends of the power emitted in free space with the one in presence of the HM one can explain the strong increase of the power ratio (

*P*

_{tot}/

*P*

_{free space}) ∝

*ω*

^{−2}in Fig. 10(b) as frequency decreases.

It is also important to provide a physical insight into the effect of the distance *h* on the power emitted by the impressed transverse dipole. In Fig. 12, we provide the plots of *P*_{tot}/*P*_{free space} and *P*_{down}/*P*_{up} at 2 THz versus the dipole’s distance *h* for the semi-infinitely periodic (*N* → ∞) graphene-dielectric multilayered structure with *μ _{c}*= 0 eV and

*μ*= 0.4 eV, obtained both via EMA (markers) and the transfer matrix method (lines). In Figs. 12(a) and 12(b) we report the power ratios

_{c}*P*

_{down}/

*P*

_{up}and

*P*

_{tot}/

*P*

_{free space}for both HMg and HMd configurations, denoting a HM with graphene (HMg) and dielectric (HMd) as topmost layer, respectively [3,49]. We notice that the responses of both HMg and HMd configurations are very similar while the HMg has a slightly larger

*P*

_{down}/

*P*

_{up}and

*P*

_{tot}/

*P*

_{free space}for smaller

*h*, in agreement with the observations in [49]. For the smallest reported distance

*h*= 0.2 μm, the emitted power (

*P*

_{tot}/

*P*

_{free space}) and

*P*

_{down}/

*P*

_{up}are maximum. However for small

*h*, EMA overestimates the reported parameters by one to two orders of magnitude for small

*h*, as demonstrated in [49], whereas for

*h*> 1 μm both methods agree well at the given frequency. Using the transfer matrix method we find the maximum ratio

*P*

_{down}/

*P*

_{up}≈ 2 × 10

^{6}when

*μ*=0 eV, and it decreases to

_{c}*P*

_{down}/

*P*

_{up}≈ 5 × 10

^{5}as a result of increasing ${Y}_{\text{down}}^{\text{TM}}$ when

*ε*′

*possesses more negative values. A similar change is also observed such that the maximum ratio*

_{t}*P*

_{tot}/

*P*

_{free space}becomes ≈ 2 × 10

^{4}when

*μ*= 0 eV, and it decreases to

_{c}*P*

_{tot}/

*P*

_{free space}≈ 4 × 10

^{2}when

*μ*= 0.4 eV. The total emitted power decreases as the distance

_{c}*h*increases, due to the stronger evanescent decay of high

*k*spectrum. However, since the distance

_{t}*h*is subwavelength, still a lot of power is able to couple into the HM. The power ratio

*P*

_{down}/

*P*

_{up}also exhibits a decrease with increasing

*h*showing the key role of the coupling of the evanescent spectrum in free space to the propagating spectrum in the HM. We finally note from the observations in Fig. 12 that the accuracy of EMA is influenced by changing the values of

*μ*, implying the effect of

_{c}*ε*on EMA’s validity.

_{t}In summary, we have shown that the power emitted by a dipolar source in the proximity of a graphene-based HM is strongly enhanced, and that it can be effectively tuned by electrostatically biasing the graphene sheets, which makes this HM a promising candidate for tunable applications in the far-infrared frequencies.

## 5. Conclusion

We have investigated a novel design of HM for far-infrared frequencies based on graphene layers. The multilayer structure has been analyzed using EMA which, based on a permittivity homogenization model, predicts the HM features at far-infrared frequencies. We have quantitatively shown the capability of tuning the composite material properties via chemical potential of graphene. We have investigated plane wave interaction with a thin film made by few graphene sheets, and showed how the transmission frequency spectrum can be tuned. We have assessed the validity of EMA for both plane wave incidence and near-field radiation from a dipole, and we have shown that under certain conditions EMA is in good agreement with the transfer matrix analysis. In the last part of the paper, we have shown that a very wide spatial spectrum emitted by an electric dipole is allowed to couple into the graphene-based HM, that would be otherwise evanescent in free space. This generates two interesting main features: (1) the power emitted by the dipole is strongly enhanced (up to several orders of magnitude) by the presence of the graphene-based HM, and (2) most of the power is directed into the HM, also for relative subwavelength HM thicknesses realized with only a few graphene sheets.These properties seem to enable the use of this tunable graphene-based HM to efficiently absorb mm-waves and terahertz frequencies, and give rise to other possible applications including super resolution lenses.

## 6. Appendix A

The transfer matrix [T_{unit}] of a unit cell composed of a graphene sheet (modeled as a lumped shunt complex admittance *Y _{s}* =

*σ*) and a dielectric layer of thickness

*d*(modeled as a transmission line) for TE

*/TM*

^{z}*waves can be written as (assuming time-harmonic fields of the form*

^{z}*e*)

^{j}^{ω}^{t}*z*axis inside the dielectric, and for TE

*and TM*

^{z}*waves: ${Z}_{d}^{\text{TE}}=\omega {\mu}_{0}/{\kappa}_{d}$ and ${Z}_{d}^{\text{TM}}={\kappa}_{d}/\left(\omega {\epsilon}_{0}{\epsilon}_{d}\right)$, respectively. We are interested in determining the Bloch wavenumber*

^{z}*k*in the

_{z}*z*-direction, that describes layer to layer propagation. Following the simple procedure in [48], the wavevector dispersion relation can be obtained from the solution of the eigenvalue problem |[T

_{unit}] −

*e*

^{−jkzd}[I]| = 0, where [I] is the identity matrix. This leads to the simple dispersion relation in Eq. (6).

## 7. Appendix B

The calculation of the admittance looking toward the −*z* direction *Y*_{HM,}* _{N}* (see Fig. 8) can be straightforwardly carried out by constructing the transfer matrix of the HM film made of

*N*unit cells [T

*], between the bottom-most material, i.e., the silicon substrate at*

_{N}*z*= −(

*Nd*+

*h*), and the surface of HM at

*z*= −

*h*. By knowing the transfer matrix of the unit cell [T

_{unit}], given in Eq. (17), one has

*Y*

_{HM,N}is evaluated using the entries of the transfer matrix [T

*] and the wave admittance inside silicon substrate,*

_{N}*Y*

_{subs}, as where ${Y}_{\text{subs}}^{\text{TE}}={\kappa}_{\text{subs}}/\left(\omega {\mu}_{0}\right)$ and ${Y}_{\text{subs}}^{\text{TM}}=\omega {\epsilon}_{0}{\epsilon}_{\text{subs}}/{\kappa}_{\text{subs}}$ are the TE

*and TM*

^{z}*wave impedances in the substrate, ${\kappa}_{\text{subs}}=\sqrt{{\epsilon}_{\text{subs}}{k}_{0}^{2}-{k}_{t}^{2}}$ is the*

^{z}*z*-directed wavenumber, and

*ε*

_{subs}=

*ε*

_{Si}is the relative permittivity of silicon. When we consider the semi-infinite case,

*N*→ ∞, the admittance

*Y*

_{HM,N}becomes the Bloch admittance

*Y*

_{Bloch}of the periodic multilayer evaluated using the unit cell’s transfer matrix entries in Eq. (17) as

*Y*

_{Bloch}such that Re(

*Y*

_{Bloch}) > 0, representing waves that carry power in the −

*z*direction.

## Acknowledgment

We thank the anonymous reviewers for their useful comments. We also would like to thank Salvatore Campione (University of California, Irvine) for fruitful discussions.

## References and links

**1. **L. Felsen and N. Marcuvitz, *Radiation and Scattering of Waves* (Prentice-Hall, NJ, 1973).

**2. **D. R. Smith and D. Schurig, “Electromagnetic wave propagation in media with indefinite permittivity and permeability tensors,” Phys. Rev. Lett. **90**, 077405 (2003) [CrossRef] [PubMed] .

**3. **C. Guclu, S. Campione, and F. Capolino, “Hyperbolic metamaterial as super absorber for scattered fields generated at its surface,” Phys. Rev. B **86**, 205130 (2012) [CrossRef] .

**4. **Y. Guo, W. Newman, C. Cortes, and Z. Jacob, “Applications of hyperbolic metamaterial substrates,” Adv. Opto-Electron. **2012**, 452502 (2012).

**5. **Y. Guo, C. L. Cortes, S. Molesky, and Z. Jacob, “Broadband super-planckian thermal emission from hyperbolic metamaterials,” Appl. Phys. Lett. **101**, 131106 (2012) [CrossRef] .

**6. **Z. Jacob, I. I. Smolyaninov, and E. E. Narimanov, “Broadband purcell effect: Radiative decay engineering with metamaterials,” Appl. Phys. Lett. **100**, 181105 (2012) [CrossRef] .

**7. **I. Smolyaninov and E. Narimanov, “Metric signature transitions in optical metamaterials,” Phys. Rev. Lett. **105**, 67402 (2010) [CrossRef] .

**8. **O. Kidwai, S. V. Zhukovsky, and J. E. Sipe, “Dipole radiation near hyperbolic metamaterials: applicability of effective-medium approximation,” Opt. Lett. **36**, 2530–2532 (2011) [CrossRef] [PubMed] .

**9. **T. U. Tumkur, J. K. Kitur, B. Chu, L. Gu, V. A. Podolskiy, E. E. Narimanov, and M. A. Noginov, “Control of reflectance and transmittance in scattering and curvilinear hyperbolic metamaterials,” Appl. Phys. Lett. **101**, 091105 (2012) [CrossRef] .

**10. **J. Pendry and S. Ramakrishna, “Refining the perfect lens,” Physica B: Condensed Matter **338**, 329 – 332 (2003) [CrossRef] .

**11. **K. J. Webb and M. Yang, “Subwavelength imaging with a multilayer silver film structure,” Opt. Lett. **31**, 2130–2132 (2006) [CrossRef] [PubMed] .

**12. **T. Tumkur, L. Gu, J. Kitur, E. Narimanov, and M. Noginov, “Control of absorption with hyperbolic metamaterials,” Appl. Phys. Lett. **100**, 161103–161103 (2012) [CrossRef] .

**13. **A. N. Poddubny, P. A. Belov, P. Ginzburg, A. V. Zayats, and Y. S. Kivshar, “Microscopic model of purcell enhancement in hyperbolic metamaterials,” Phys. Rev. B **86**, 035148 (2012) [CrossRef] .

**14. **C. L. Cortes, W. Newman, S. Molesky, and Z. Jacob, “Quantum nanophotonics using hyperbolic metamaterials,” J. Opt. **14**, 063001 (2012) [CrossRef] .

**15. **G. V. Naik, J. Liu, A. V. Kildishev, V. M. Shalaev, and A. Boltasseva, “Demonstration of al:zno as a plasmonic component for near-infrared metamaterials,” PNAS **109**, 8834–8838 (2012) [CrossRef] [PubMed] .

**16. **J. Kim, V. Drachev, Z. Jacob, G. Naik, A. Boltasseva, E. Narimanov, and V. Shalaev, “Improving the radiative decay rate for dye molecules with hyperbolic metamaterials,” Opt. Express **20**, 8100–8116 (2012) [CrossRef] [PubMed] .

**17. **G. Naik and A. Boltasseva, “Semiconductors for plasmonics and metamaterials,” Phys. Status Solidi Rapid Res. Lett. **4**, 295–297 (2010) [CrossRef] .

**18. **C. Rizza, A. Ciattoni, E. Spinozzi, and L. Columbo, “Terahertz active spatial filtering through optically tunable hyperbolic metamaterials,” Opt. Lett. **37**, 3345–3347 (2012) [CrossRef] .

**19. **A. Vakil and N. Engheta, “One-atom-thick reflectors for surface plasmon polariton surface waves on graphene,” Opt. Comm. **285**, 3428 – 3430 (2012) [CrossRef] .

**20. **F. Rana, “Graphene terahertz plasmon oscillators,” IEEE Trans. Nanotechnol. **7**, 91–99 (2008) [CrossRef] .

**21. **G. W. Hanson, A. B. Yakovlev, and A. Mafi, “Excitation of discrete and continuous spectrum for a surface conductivity model of graphene,” J. Appl. Phys. **110**, 114305 (2011) [CrossRef] .

**22. **M. Tamagnone, J. Gomez-Diaz, J. Mosig, and J. Perruisseau-Carrier, “Analysis and design of terahertz antennas based on plasmonic resonant graphene sheets,” J. Appl. Phys. **112**, 114915 (2012) [CrossRef] .

**23. **B. Wang, X. Zhang, F. J. Garcia-Vidal, X. Yuan, and J. Teng, “Strong coupling of surface plasmon polaritons in monolayer graphene sheet arrays,” Phys. Rev. Lett. **109**, 073901 (2012) [CrossRef] [PubMed] .

**24. **I. V. Iorsh, I. S. Mukhin, I. V. Shadrivov, P. A. Belov, and Y. S. Kivshar, “Hyperbolic metamaterials based on multilayer graphene structures,” Phys. Rev. B **87**, 075416 (2013) [CrossRef] .

**25. **A. Andryieuski, A. V. Lavrinenko, and D. N. Chigrin, “Graphene hyperlens for terahertz radiation,” Phys. Rev. B **86**, 121108 (2012) [CrossRef] .

**26. **V. P. Gusynin and S. G. Sharapov, “Unconventional integer quantum hall effect in graphene,” Phys. Rev. Lett. **95**, 146801 (2005) [CrossRef] [PubMed] .

**27. **A. B. Kuzmenko, E. van Heumen, F. Carbone, and D. van der Marel, “Universal optical conductance of graphite,” Phys. Rev. Lett. **100**, 117401 (2008) [CrossRef] [PubMed] .

**28. **L. Gerhard, E. Moyen, T. Balashov, I. Ozerov, M. Portail, H. Sahaf, L. Masson, W. Wulfhekel, and M. Hanbucken, “A graphene electron lens,” Appl. Phys. Lett. **100**, 153106 (2012) [CrossRef] .

**29. **A. Vakil and N. Engheta, “Transformation optics using graphene,” Science **332**, 1291–1294 (2011) [CrossRef] [PubMed] .

**30. **C. Chen, S. Rosenblatt, K. Bolotin, W. Kalb, P. Kim, I. Kymissis, H. Stormer, T. Heinz, and J. Hone, “Performance of monolayer graphene nanomechanical resonators with electrical readout,” Nature Nanotech. **4**, 861–867 (2009) [CrossRef] .

**31. **X. Wang, L. Zhi, and K. Mullen, “Transparent, conductive graphene electrodes for dye-sensitized solar cells,” Nano Lett. **8**, 323–327 (2008) [CrossRef] .

**32. **C. S. R. Kaipa, G. W. P. Y. R. Yakovlev, Alexander Hanson, M. F. Medina, and F., “Enhanced transmission with a graphene-dielectric microstructure at low-terahertz frequencies,” Phys. Rev. B **85**, 245407 (2012) [CrossRef] .

**33. **S. Thongrattanasiri, F. H. L. Koppens, and F. J. Garcia de Abajo, “Complete optical absorption in periodically patterned graphene,” Phys. Rev. Lett. **108**, 047401 (2012) [CrossRef] [PubMed] .

**34. **A. Fallahi and J. Perruisseau-Carrier, “Design of tunable biperiodic graphene metasurfaces,” Phys. Rev. B **86**, 195408 (2012) [CrossRef] .

**35. **D. Sounas and C. Caloz, “Gyrotropy and nonreciprocity of graphene for microwave applications,” IEEE Trans. Microw. Theory Techn. **60**, 901 –914 (2012) [CrossRef] .

**36. **D. L. Sounas and C. Caloz, “Electromagnetic nonreciprocity and gyrotropy of graphene,” Appl. Phys. Lett. **98**, 021911 (2011) [CrossRef] .

**37. **G. Lovat, “Equivalent circuit for electromagnetic interaction and transmission through graphene sheets,” IEEE Trans. Electromagn. Compat. **54**, 101 –109 (2012) [CrossRef] .

**38. **J. Sun, J. Zhou, B. Li, and F. Kang, “Indefinite permittivity and negative refraction in natural material: Graphite,” Appl. Phys. Lett. **98**, 101901 (2011) [CrossRef] .

**39. **C. S. Kaipa, A. B. Yakovlev, F. Medina, F. Mesa, C. Butler, and A. P. Hibbins, “Circuit modeling of the transmissivity of stacked two-dimensional metallic meshes,” Opt. Express **18**, 13309–13320 (2010) [CrossRef] [PubMed] .

**40. **Y. R. Padooru, A. B. Yakovlev, C. S. Kaipa, F. Medina, and F. Mesa, “Circuit modeling of multiband high-impedance surface absorbers in the microwave regime,” Phys. Rev. B **84**, 035108 (2011) [CrossRef] .

**41. **C. S. R. Kaipa, A. B. Yakovlev, F. Medina, and F. Mesa, “Transmission through stacked 2d periodic distributions of square conducting patches,” J. Appl. Phys. **112**, 033101 (2012) [CrossRef] .

**42. **X. Ni, G. Naik, A. Kildishev, Y. Barnakov, A. Boltasseva, and V. Shalaev, “Effect of metallic and hyperbolic metamaterial surfaces on electric and magnetic dipole emission transitions,” Appl. Phys. B **103**, 553–558 (2011) [CrossRef] .

**43. **K. Novoselov, A. Geim, S. Morozov, D. Jiang, Y. Zhang, S. Dubonos, I. Grigorieva, and A. Firsov, “Electric field effect in atomically thin carbon films,” Science **306**, 666–669 (2004) [CrossRef] [PubMed] .

**44. **V. P. Gusynin, S. G. Sharapov, and J. P. Carbotte, “Sum rules for the optical and hall conductivity in graphene,” Phys. Rev. B **75**, 165407 (2007) [CrossRef] .

**45. **G. W. Hanson, “Dyadic green’s functions and guided surface waves for a surface conductivity model of graphene,” J. Appl. Phys. **103**, 064302 (2008) [CrossRef] .

**46. **R. A. Jishi, M. S. Dresselhaus, and G. Dresselhaus, “Electron-phonon coupling and the electrical conductivity of fullerene nanotubules,” Phys. Rev. B **48**, 11385–11389 (1993) [CrossRef] .

**47. **Y.-W. Tan, Y. Zhang, K. Bolotin, Y. Zhao, S. Adam, E. H. Hwang, S. Das Sarma, H. L. Stormer, and P. Kim, “Measurement of scattering rate and minimum conductivity in graphene,” Phys. Rev. Lett. **99**, 246803 (2007) [CrossRef] .

**48. **D. Pozar, *Microwave engineering* (John Wiley & Sons, 2009).

**49. **O. Kidwai, S. V. Zhukovsky, and J. E. Sipe, “Effective-medium approach to planar multilayer hyperbolic meta-materials: Strengths and limitations,” Phys. Rev. A **85**, 053842 (2012) [CrossRef] .

**50. **H. N. S. Krishnamoorthy, Z. Jacob, E. Narimanov, I. Kretzschmar, and V. M. Menon, “Topological transitions in metamaterials,” Science **336**, 205–209 (2012) [CrossRef] [PubMed] .

**51. **S. Campione, S. Steshenko, M. Albani, and F. Capolino, “Complex modes and effective refractive index in 3d periodic arrays of plasmonic nanospheres,” Opt. Express **19**, 26027–26043 (2011) [CrossRef] .

**52. **F. Capolino and M. Albani, “Truncation effects in a semi-infinite periodic array of thin strips: A discrete wiener-hopf formulation,” Radio Sci. **44**, RS2S91 (2009) [CrossRef] .

**53. **P.-Y. Chen and A. Alu, “Atomically thin surface cloak using graphene monolayers,” ACS Nano **5**, 5855–5863 (2011) [CrossRef] [PubMed] .

**54. **Z. Jacob, J. Kim, G. Naik, A. Boltasseva, E. Narimanov, and V. Shalaev, “Engineering photonic density of states using metamaterials,” Appl. Phys. B **100**, 215–218 (2010) [CrossRef] .