## Abstract

We analyze and model the nonlocal response of ultrathin hyperbolic metasurfaces (HMTSs) by applying an effective medium approach. We show that the intrinsic spatial dispersion in the materials employed to realize the metasurfaces imposes a wavenumber cutoff on the hyperbolic isofrequency contour, inversely proportional to the Fermi velocity, and we compare it with the cutoff arising from the structure granularity. In the particular case of HTMSs implemented by an array of graphene nanostrips, we find that graphene nonlocality can become the dominant mechanism that closes the hyperbolic contour – imposing a wavenumber cutoff at around 300${k}_{0}$ – in realistic configurations with periodicity $L<\pi /(300{k}_{0})$, thus providing a practical design rule to implement HMTSs at THz and infrared frequencies. In contrast, more common plasmonic materials, such as noble metals, operate at much higher frequencies, and therefore their intrinsic nonlocal response is mainly relevant in hyperbolic metasurfaces and metamaterials with periodicity below a few nm, being very weak in practical scenarios. In addition, we investigate how spatial dispersion affects the spontaneous emission rate of emitters located close to HMTSs. Our results establish an upper bound set by nonlocality to the maximum field confinement and light-matter interactions achievable in practical HMTSs, and may find application in the practical development of hyperlenses, sensors and on-chip networks.

© 2015 Optical Society of America

## 1. Introduction

Hyperbolic metasurfaces (HMTSs) are extremely anisotropic 2D structures whose transverse surface conductivity changes the sign of its imaginary part as a function of the polarization of the in-plane electric field [1], thus behaving as a dielectric along one direction and as a metal along the orthogonal one within the sheet. The ultrathin planar nature of HTMSs solves some of the major practical challenges of bulk hyperbolic metamaterials (HMTMSs) [2–4], thanks to a large resilience to volumetric losses, a relatively simple fabrication using standard lithographic and etching techniques [5], and easy external access to the propagating energy, while being fully compatible with integrated circuits and optoelectronic components. In contrast to HMTMs, HMTSs support the propagation of low-loss and extremely confined surface plasmon polaritons (SPPs) with evanescent fields in the direction perpendicular to the surface, providing access to a drastic enhancement of light-matter interactions on a surface [1]. Furthermore, HMTSs allow a large degree of freedom to manipulate the supported SPPs, including routing them towards specific directions within the sheet, dispersion-free propagation (canalization), and negative refraction [5]. The intriguing electromagnetic response of HMTSs arises from their hyperbolic isofrequency contour, associated with an ideally infinite local density of states (LDOS). In classical terms, LDOS and the maximum guided wave number are limited by two mechanisms, namely the presence of losses [1] and the granularity of the metasurface [1, 6], that impose a wavenumber cutoff and close the otherwise open hyperbolic contour. In the case of HMTMs, a third fundamental mechanism has recently been revealed by considering the actual nonlocal response of the metals composing the structure [7], a phenomenon that gives rise to a wavenumber cutoff inversely proportional to the Fermi velocity of electrons in the material [7].

There are several possibilities to implement HMTSs in practice. At visible frequencies, HMTSs based on single-crystalline silver nanostructures have been experimental demonstrated in [5]. Another approach relies on using anisotropic subwavelength metallic scatterers, as envisaged in [8]. An interesting alternative operating at THz and near infrared frequencies uses a densely packed array of graphene nanostrips [1]. This configuration – analyzed using an effective medium approach (EMA) in [1,6] and the Kubo formalism in [9]– allows tuning in real time the wavenumber contour topology, from closed elliptical to open hyperbolic going through the extremely anisotropic $\sigma $-near zero case [1], by simply applying a gate bias. This structure takes advantage of the inherent reconfigurable capabilities of graphene plasmonics [10,11], a field that has recently emerged as an excellent platform for strong light-matter interactions and has triggered a myriad of exciting applications [12–16]. In this context, the intrinsic spatial dispersion of graphene [17–19] is known to significantly affect the response of guided THz waveguides and components [20–22] and to blueshift the plasmonic resonances appearing in graphene nanostructures [23].

In this contribution, we investigate the electromagnetic nonlocal response of hyperbolic metasurfaces, focusing on the case of a densely-packed array of graphene strips at THz and near infrared frequencies. We have chosen this configuration because the intrinsic nonlocal effects of graphene may be stronger than those usually found in noble metals, and they can be strong enough to become the fundamental mechanism that closes the hyperbolic dispersion relation of the structure. However, we also consider implementations based on other 2D materials, assuming that their nonlocal conductivity is available. Our main goal is to isolate and independently investigate nonlocalities arising due to (i) granularity and (ii) the intrinsic spatially-dispersive response of the constituent materials. The former case is analyzed using a full-wave mode-matching approach [24], where graphene is characterized using a local model. As expected, the lattice periodicity *L* closes the hyperbolic isofrequency contour of the structure by introducing a cutoff wavenumber at $\approx \pi /L.$ In the latter scenario, we isolate the influence of graphene nonlocality in the response of the entire HMTSs by applying an effective medium approach that homogenizes the 2D structure in the $L\to 0$ limit, thus removing the spatial dispersion introduced by the lattice. This allows us to demonstrate that the inherent nonlocal response of graphene imposes a cutoff to the HMTSs wavenumbers inversely proportional to the Fermi velocity of the material, as occurs in HMTMs [7]. We do stress that graphene possesses a relatively high Fermi velocity (${v}_{F}\approx {10}^{6}$ m/s) able to significantly modify the propagation of surface waves with ${k}_{\rho}>300{k}_{0}$. Consequently, graphene nonlocality becomes the dominant mechanism that controls SPPs propagation in hyperbolic metasurfaces with periodicity $L<\pi /300{k}_{0}$, thus setting a simple rule for the design of practical HMTSs. Taking into account the operation frequency of these structures, nonlocal effects play a key role in the response of realistic HMTSs with periodicity ranging between 500 to 20 nm. This result is in clear contrast to the case of conventional optical HMTMs, where intrinsic spatial dispersion effects of the materials may be relevant only for subwavelength composites – with periodicity below a few nm –very challenging to produce even with advanced nanofabrication techniques. Therefore, they are usually of limited practical impact. We conclude our study showing that nonlocality affects the spontaneous emission rate (SER) of emitters located very close to HMTSs, imposing a fundamental limit associated to the large wavevector cutoff. Throughout the paper, full-wave simulations based on the commercial software COMSOL Multiphysics are employed to validate our analytical study.

## 2. Homogenization of nonlocal HMTSs

In this section, we first describe the hyperbolic metasurface under study – a densely packed array of 2D strips – and derive an EMA able to homogenize the structure in the general case of ribbons characterized by a fully-populated conductivity tensor. Then, we introduce and discuss spatially dispersive models employed to characterize graphene and 2D materials. Next, we study the effective conductivity of homogenized nonlocal HTMSs. Finally, we derive the dispersion relation of theses structures and provide brief guidelines to its numerical solution.

#### 2.1 EMA of periodic 2D strips defined by a fully-populated conductivity tensor

Let us consider an array of densely packed 2D strips, as illustrated in Fig. 1, where each strip is described by a fully populated conductivity tensor

*L*is the periodicity of the unit-cell, and

*G*is the separation distance between two adjacent strips. This description allows transforming the original structure into an equivalent one composed of strips with different conductivities, as shown in Fig. 1.

The main goal of this subsection is to derive an EMA to characterize the effective conductivity tensor of the metasurface in the limit of infinitesimally small periodicity, i.e., when $L\to 0$. To this purpose, we extend the EMA introduced in [1,26]. Specifically, applying boundary conditions to the equivalent unit-cell (see inset of Fig. 1) we note that (i) the electric field ${E}_{y}$ must be continuous at the strip edges, and (ii) there is no specific restriction that prevents the continuity of the surface current across the ribbons. Taking advantage of the subwavelength nature of the structure, and given the boundary conditions on the parallel strips (continuous normal current distribution and continuous tangential electric field at the strip edges), we enforce a constant current ${J}_{x}$ and electric field ${E}_{y}$ across theentire unit cell. Then, applying Ohm’s law, one finds

#### 2.2 Intrinsic nonlocal response of 2D materials

We analyze here the influence of nonlocal effects in the electromagnetic response of 2D bare materials. We focus on graphene at THz and near infrared frequencies, a material that has recently led to the development of exciting HTMSs. Then, we briefly discuss the possible nonlocal response of usual metals assuming ultrathin (2D) configurations. Throughout this work, we model spatially-dispersive graphene using the Bhatnagar-Gross-Krook (BGK) approach derived in [19]. This model is based on the Boltzmann transport equation, taking into account intraband contributions in graphene, and it is accurate up to tens of THz when the spatial variations of the fields are lower than the de Broglie wavelength of the particles [19] (i.e. ${k}_{\rho}<2{k}_{F}$, where ${k}_{F}$ is the Fermi wavenumber). According to this model, graphene is modelled by a fully-populated conductivity tensor $\overline{\overline{{\sigma}^{g}}}$ defined by

*isotropic*in the

*xy*plane [19].

Although we are mainly interested in graphene-based HMTSs in this work, it might be instructive to analyze the possible nonlocal conductivity of other 2D metals. For instance, assuming for a moment that the hydrodynamic Drude model within the Thomas-Fermi approximation [7, 32] holds as we thin down a metal and that we operate above the plasma frequency, its effective in-plane conductivity can be expressed in polar coordinates as

*t*’ is the metal thickness, and

#### 2.3 Nonlocal effective parameters of homogenized HMTSs

The effective conductivity tensor $\overline{\overline{{\sigma}^{eff}}}\left({k}_{x},{k}_{y}\right)$ of nonlocal HMTSs composed of densely packed 2D strips can easily be derived by combining the EMA derived above with the spatially-dispersive model of the materials that compose the structure. In the general case, it can be expressed as

*xx*component of the effective conductivity tensor is dominated by the near-field coupling between the strips, thus leading to an almost constant dielectric (capacitive) behavior for all wavenumbers. However, spatial dispersion may significantly affect the response of the other conductivity components. In the case of low wavenumbers, ${\sigma}_{yy}^{eff}$ keeps its local metallic nature and the cross-terms are negligible, thus leading to the usual description of hyperbolic metasurfaces [1].

This scenario is different as the wavenumbers approach $\left|{k}_{\rho}/{k}_{0}\right|\approx c/{v}_{F}\approx 300$, where the transverse part of nonlocal graphene conductivity dominates the response of the material. As a consequence, ${\sigma}_{yy}^{eff}$ changes its nature from metallic to dielectric for large ${k}_{y}$ wavenumbers, while the cross-terms also contributes to the response. Intuitively, this behavior can be better understood in polar coordinates, since the nonlocal longitudinal conductivity is strongly reduced –and changes the nature of its response from metallic to dielectric–preventing energy propagation along the strips, whereas the dominant transverse nonlocal component rotates the energy flow towards the ribbons edges. We emphasize again that the behavior described above is not by any means restricted to the case of graphene, but it also applies to other HMTSs composed of nonlocal 2D materials, such as those described by the hydrodynamic Drude model within the Thomas-Fermi approximation.

#### 2.4 Dispersion relation of homogenized nonlocal HMTSs

The dispersion relation of homogeneous ultrathin metasurfaces defined by an effective conductivity tensor $\overline{\overline{{\sigma}^{eff}}}.$ can be expressed as [33]

In the local case, Eq. (21) can be solved analytically by fixing the direction of propagation of the supported SPPs and then physically rotating the metasurface, as described in [6]. However, in the nonlocal case this approach cannot be applied due to the intrinsic spatial dispersion of $\overline{\overline{{\sigma}^{eff}}}$, and one has to resort to numerical methods. In this work, we solve Eq. (21) by applying the Newton–Raphson approach [34]. In order to avoid finding local minima in the complex plane, and to speed up the convergence of the algorithm, it is in general most efficient to first find the solution in the *y*-direction, where the hyperbolic plasmons are less confined and the local dispersion relation provides a good starting point, and then progressively track the solution in other directions until cutoff.

## 3. Influence of nonlocalities on the electromagnetic response of HTMSs

In this Section, we discuss the influence of different nonlocalities on the isofrequency contours of HTMSs and on the SER of emitters located in their vicinity, and we determine under which conditions they may dominate the electromagnetic response of the entire structure. To this purpose, we apply the local and nonlocal versions of the EMA derived in Section II, comparing the predicted dispersion relations with results computed by two independent full-wave methods. First, to isolate the influence of spatial dispersion associated with the macroscopic periodicity, we utilize a mode matching approach that considers a local, scalar ${\sigma}^{g}$, and is able to accurately predict the closing of the hyperbolic contour due to lattice granularity [1, 6, 24]. Then, in order to include the intrinsic nonlocal response of graphene, we use the commercial FEM solver COMSOL Multiphysics to find the complex eigenfrequencies of the unit cell for arbitrary macroscopic wavevectors. Here, graphene is modelled as a surface current with the BGK model described above, thus taking into account both types of spatial dispersion. However, we stress that our COMSOL simulations are only strictly valid in the $L\to 0$ limit, when the spatial variation of the fields within the unit cell is only associated to the macroscopic propagation along the metasurface, and provide a very good approximation otherwise. A completely rigorous numerical solution would require the development of an additional 2D boundary condition, as similarly done for nonlocal bulk media [7, 30, 31]. We do remark that this approach is justified here, since we are mainly interested in finding the inherent limitations of nonlocal HMTSs even when structure granularity can be neglected. Moreover, in the case of larger unit cells, the hyperbolic contour will be closed before the nonlocal conductivity becomes relevant, as it will be shown below.Figure 5 shows the isofrequency contour of an HMTS with a fixed graphene filling factor of 0.5, operating at $f=2.5$ THz, with ${\mu}_{c}=0.02$ eV and $\tau =0.3$ ps, for several values of the lattice period *L*. This approach allows us to easily identify the different physical mechanisms responsible for closing the otherwise open hyperbolic isofrequency contour. In case of extremely subwavelength unit cells [Fig. 5(a), $L=30$ nm], the lattice periodicity imposes a very large wavenumber cutoff, and therefore the approaches based on local material predict a very similar response. In this scenario, the intrinsic nonlocal response of graphene enforces a wavenumber cutoff at around $\left|{k}_{\rho}\right|\approx (c/{v}_{F}){k}_{0}$, fully consistent with the arguments given in Section II. Importantly, the supported SPPs abruptly disappear instead of following the expected closing towards the ${k}_{x}$ axis. This behavior arises because at large wavenumbers the transverse component of nonlocal graphene conductivity dominates, thus forcing energy towards the ribbons edges and effectively rotating the SPPs angle of propagation within the sheet. In addition, we do note that HMTSs only supports propagation through a well-defined set of directions, as described in [8]. The combination of these two factors fully explains the mode disappearance in the nonlocal isofrequency contour of Fig. 5(a). We stress the excellent agreement found between our nonlocal EMA approach and full-wave simulations using COMSOL. If we consider now $L=100$ nm [Fig. 5(b)], the graphene spatial dispersion remains the dominant mechanism, but the structure granularity becomes non-negligible near cutoff, as shown by the slight deviation of the numerical results from the predicted EMA dispersion. In this case, full-wave simulations using local conductivity are still very inaccurate compared to the simple EMA developed here. A complementary scenario is shown in Fig. 5(c) ($L=300$ nm), where granularity is the main mechanism responsible of the closing of the contour, but graphene spatial dispersion further deforms it towards the *x* direction. In this case, both nonlocal phenomena should be taken into account. Finally, Fig. 5(d) considers $L=700$ nm. Now, graphene nonlocality is negligible, since the hyperbola is closed at low wavenumbers due to the large lattice period.

Our analysis of nonlocal graphene-based HMTSs provides important design guidelines for practical devices. Specifically, we have demonstrated that reducing the size of the unit cell in order to increase the LDOS is only beneficial until $L\approx {v}_{F}\pi /\left(c{k}_{0}\right)$, since the inherent nonlocal response of electrons becomes dominant at this point. Therefore, nonlocality imposes a practical limit to the HMTSs periodicity that leads to realizable structures at THz and near infrared frequencies - fully within the capabilities of current nanofabrication technology [36] - able to reach the maximum possible wave confinement and LDOS.

Once we have validated the accuracy of our proposed nonlocal EMA, we apply it to illustrate the response of graphene-based HMTSs versus some physical parameters of the structure. First, we investigate in Fig. 6(a) the influence of the strip width *W*, for a fixed period *L*. Interestingly, the slightly open branches associated to a hyperbolic dispersion with reduced ${\sigma}_{yy}$ (small *W/L*, red lines) obtained using a local EMA are compensated by graphene nonlocality, thus leading to an almost perfect canalization. For other strip widths, nonlocal effects shows a relatively minor influence on the hyperbolic contour until cutoff. Figure 6(b) illustrates the robustness of the hyperbolic response of spatially dispersive HMTSs versus graphene loss, similarly to previous studies that considered local conductivity [1]. We note, however, that the effective propagation length of the supported hyperbolic plasmons is strongly dependent on graphene therelaxation time $\tau $. In this sense, we have observed that the figure of merit ratio $Re\left[{k}_{\rho}\right]/Im\left[{k}_{\rho}\right]$ is slightly lower than the case of bare graphene [1, 6]. Finally, Figs. 6(c)-6(d) showcase the reconfiguration capabilities of HMTS through the control of graphene chemical potential ${\mu}_{c}$, easily achieved through electrostatic biasing of the strips [11]. We consider two scenarios: Fig. 6(c) corresponds to an HMTS designed to operate far from potential topological transitions, with a hyperbolic dispersion that is widely tunable. Spatial dispersion limits the maximum confinement in all cases, most noticeably for ${\mu}_{c}=0.01$ eV due to the larger minimum wavenumber supported. On the other hand, Fig. 6(d) depicts the isofrequency contour of a HMTS operating in the ${\sigma}_{xx}$ near-zero regime, where small modifications of the chemical potential may result in topological transitions from elliptical to hyperbolic. Interestingly, this study demonstrates *nonlocality-induced topological transitions* [35]. For ${\mu}_{c}=0.01\text{}\text{eV}$ both local and nonlocal EMA models predict a closed contour. However, spatial dispersion causes the transition to hyperbolic topology to occur for different values of ${\mu}_{c}$, as shown by the plots for ${\mu}_{c}=0.075$ eV and ${\mu}_{c}=0.2$ eV. Consequently, while the local model predicts a closed contour, nonlocality forces a topological transition that leads to open hyperbolic dispersion. This behavior further demonstrates the importance of correctly modelling nonlocality in the design of HTMSs, since it shows that the local description of the materials may predict an incorrect topology.

Enhancement of light-matter interactions is one of the most attractive features of HMTSs [1, 6], so it is indeed of interest to study how spatial dispersion affects the SER of emitters located in their vicinity. To this end, we calculate the SER of a *z*-oriented dipole placed above a graphene-based HMTSs with $L=50$ nm, $W=25$ nm, ${\mu}_{c}$ = 0.02 eV, $\tau =0.3$ ps, and$f=2.5$ THz. Results are obtained following the approach recently introduced in [1, 6], modified here to include graphene nonlocal response. Figures 7(a)-7(b) show the emitter SER versus the ribbon width and the distance *d* between the source and the surface, applying the local and nonlocal EMAs. In order to reveal the influence of nonlocality to the SER, we focus here on emitters closely located to HTMSs, thus preventing free space from filtering waves with high spatial frequency. For local graphene, very high SER is predicted close to the canalization regime – implemented by very narrow ribbons – due to the extremely high light confinement enabled by this topology [1,6], up to thousands of times larger than ${k}_{0}$. Essentially, this design maximizes the range of ${k}_{x}$ components radiated by the emitter that can be guided by the metasurface. Since we are considering small distances between the emitter and the HMTS, these components are not fully filtered by free space and a large portion of their power is coupled to the HMTS. In the case of nonlocal graphene, however, components with ${k}_{\rho}>300{k}_{0}$ are forbidden, leading to a reduced SER. Contrary to the local case, Fig. 7(b) shows that for small distances, spatial dispersion causes the SER to increase with *W*, reaching an optimum at 41 nm. We have verified that this *W* corresponds precisely to the design that supports the widest range of wavenumbers, consistent with the conclusions in the nondispersive scenario. Note that for larger separations between emitter and HMTS both approaches leads to very similar results, due to the filtering of evanescent waves by free-space. Interestingly, these results also show that spatial dispersion does not necessarily cause a reduction of the SER for every HMTS design, although it does impose upper limits associated to the high wavevector cutoff.

In order to further investigate this limit, we investigate in Fig. 8(a) the SER of a *z*-oriented dipole versus its distance to a lossless homogenized metasurface. We have removed in this example the presence of losses, since the standard definition of SER [39] does not hold when the emitter in embedded in lossy media [40, 41]. As expected [7], this leads to a significant decrease of the emitter SER compared to the examples shown in Figs. 7(a)-7(b). In the local case, neglecting the influence of the granularity, the emitter SER diverges as the emitter gets closer to the HMTSs, ideally providing an infinite SER when the dipole is placed *exactly* on the sheet. The picture is drastically different when the nonlocal response of graphene is taken into account, since it removes the singularity and provides a finite and realistic response. This example clearly demonstrates how nonlocality imposes a fundamental limit to the response of hyperbolic metasurfaces. The features of this limit mainly depend on the Fermi velocity of theelectrons within the nonlocal 2D material, and increases as the strip width or chemical potential is reduced [see Fig. 8(b)]. As expected, emitters located on structures that support surface waves with elliptic dispersion provide SER orders of magnitude lower than in the hyperbolic case.

Finally, we would like to emphasize that the qualitative behavior of nonlocal HMTSs discussed in this section and the fundamental limits that they impose on the LDOS and wave confinement arise by the finite Fermi velocity of electrons, and therefore they are not restricted to graphene-based nanostructures, but they hold for any HMTS.

## 4. Conclusions

In conclusion, we have investigated the nonlocal response of HMTSs in detail. To this purpose, we have derived an accurate EMA to homogenize nonlocal 2D metasurfaces assuming an infinitesimally small periodicity $\left(L\to 0\right)$. This approach has allowed us to remove the spatial dispersion effects introduced by the lattice in order to focus on the influence of the material nonlocality in the response of HMTSs. We do note that our approach constitutes a first approximation to the rigorous analysis of nonlocal HTMSs. More advanced models should take into account: (i) the influence of the lattice in the homogenization process [37], and (ii) an additional boundary condition able to impose the continuity of the tangential surface current in 2D media with intrinsic spatial dispersion [7, 30, 31]. Applying our EMA, we have demonstrated that material nonlocality imposes a wavenumber cutoff to the open isofrequency contour of HMTSs that is inversely proportional to the electron Fermi velocity, a phenomenon also found in bulk HMTMs [7]. We have then compared this wavenumber cutoff to the usual one arising from the metasurface granularity, delimiting the conditions under which each type of nonlocality dominates the HMTSs response. Finally, we have studied the influence of spatial dispersion in the SER of emitters located very close to HMTSs.

Our results are particularly relevant in case of graphene-based HMTSs operating at THz and near infrared frequencies. There, the intrinsic spatial dispersion of graphene becomes the dominant mechanism to close the hyperbolic contour – at around ${k}_{\rho}\approx 300{k}_{0}$ – in practical implementations with periodicities $L<\pi /(300{k}_{0})$, i.e., in the range of 500-20 nm. This result clearly illustrates the tremendous importance of nonlocal effects in the correct analysis and design of realistic HMTSs. We stress that similar graphene-based nanostructures have already been designed and fabricated for other purposes [36, 38]. These configurations may therefore provide an excellent scenario to indirectly measure graphene nonlocality. In the case of conventional HMTSs and HMTMs at optical wavelengths, we have found that the intrinsic nonlocal effects of metals are only relevant in configurations with periodicity below a few nanometers. Since these devices are usually very challenging to produce, even with modern nanofabrication techniques, this type of spatial dispersion is usually expected to be very weak in practice. Even though the intrinsic nonlocal effects of materials impose upper bounds to the field confinement and LDOS achievable in HMTSs, there is still plenty of room for exciting applications in realistic structures operating from THz to optics, including planar hyperlenses, imaging systems, sensors exploiting the extremely large (but finite) light-matter interactions, and on-chip networks.

## Acknowledgments

The authors wish to thank Dr. D. Sounas (The University of Texas of Austin, US) for fruitful discussions. This work was supported by the Air Force Office of Scientific Research (AFOSR) with grant No. FA9550-13-1-0204, the Welch foundation with grant No. F-1802, and the National Science Foundation with grant No. ECCS-1406235.

## References and links

**1. **J. S. Gomez-Diaz, M. Tymchenko, and A. Alù, “Hyperbolic Plasmons and Topological Transitions Over Uniaxial Metasurfaces,” Phys. Rev. Lett. **114**(23), 233901 (2015). [CrossRef] [PubMed]

**2. **A. Poddubny, I. Iorsh, P. Belov, and Y. Kivshar, “Hyperbolic metamaterials,” Nat. Photonics **7**(12), 948–957 (2013). [CrossRef]

**3. **V. P. Drachev, V. A. Podolskiy, and A. V. Kildishev, “Hyperbolic metamaterials: new physics behind a classical problem,” Opt. Express **21**(12), 15048–15064 (2013). [CrossRef] [PubMed]

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

**5. **A. A. High, R. C. Devlin, A. Dibos, M. Polking, D. S. Wild, J. Perczel, N. P. de Leon, M. D. Lukin, and H. Park, “Visible-frequency hyperbolic metasurface,” Nature **522**(7555), 192–196 (2015). [CrossRef] [PubMed]

**6. **J. S. Gomez-Diaz, M. Tymchenko, and A. Alù, “Hyperbolic metasurfaces: surface plasmons, light-matter interactions, and physical implementation using graphene strips,” Opt. Mater. Express **5**(10), 2313–2329 (2015). [CrossRef]

**7. **W. Yan, M. Wubs, and N. A. Mortensen, “Hyperbolic metamaterials: Nonlocal response regularizes broadband supersingularity,” Phys. Rev. B **86**(20), 205429 (2012). [CrossRef]

**8. **O. Y. Yermakov, A. I. Ovcharenko, M. Song, A. A. Bogdanov, I. V. Iorsh, and Yu. S. Kivshar, “Hybrid waves localized at hyperbolic metasurfaces,” Phys. Rev. B **91**(23), 235423 (2015). [CrossRef]

**9. **I. Trushkov and I. Iorsh, “Two-dimensional hyperbolic medium for electrons and photons based on the array of tunnel-coupled graphene nanoribbons,” Phys. Rev. B **92**(4), 045305 (2015). [CrossRef]

**10. **F. H. L. Koppens, D. E. Chang, and F. J. García de Abajo, “Graphene Plasmonics: A Platform for Strong Light-Matter Interactions,” Nano Lett. **11**(8), 3370–3377 (2011). [CrossRef] [PubMed]

**11. **F. J. García de Abajo, “Graphene Plasmonics: Challenges and Opportunities,” ACS Photonics **1**(3), 135–152 (2014). [CrossRef]

**12. **J. S. Gómez-Díaz and J. Perruisseau-Carrier, “Graphene-based plasmonic switches at near infrared frequencies,” Opt. Express **21**(13), 15490–15504 (2013). [CrossRef] [PubMed]

**13. **J. S. Gómez-Díaz, M. Esquius-Morote, and J. Perruisseau-Carrier, “Plane wave excitation-detection of non-resonant plasmons along finite-width graphene strips,” Opt. Express **21**(21), 24856–24872 (2013). [CrossRef] [PubMed]

**14. **P. Chen, H. Huang, D. Akinwande, and A. Alù, “Graphene-Based Plasmonic Platform for Reconfigurable Terahertz Nanodevices,” ACS Photonics **1**(8), 647–654 (2014). [CrossRef]

**15. **D. L. Sounas, H. S. Skulason, H. V. Nguyen, A. Guermoune, M. Siaj, T. Szkopek, and C. Caloz, “Faraday rotation in magnetically biased graphene at microwave frequencies,” Appl. Phys. Lett. **102**(19), 191901 (2013). [CrossRef]

**16. **S. Thongrattanasiri, F. H. L. Koppens, and F. J. García de Abajo, “Complete Optical Absorption in Periodically Patterned Graphene,” Phys. Rev. Lett. **108**(4), 047401 (2012). [CrossRef] [PubMed]

**17. **L. A. Falkovsky and A. A. Varlamov, “Space-time dispersion of graphene conductivity,” Eur. Phys. J. B **56**(4), 281–284 (2007). [CrossRef]

**18. **M. Jablan, H. Buljan, and M. Soljacic, “Plasmonics in graphene at infrared frequencies,” Phys. Rev. B Condens. Matter **80**(24), 245435 (2009). [CrossRef]

**19. **G. Lovat, G. W. Hanson, R. Araneo, and P. Burghignoli, “Semiclassical spatially dispersive intraband conductivity tensor and quantum capacitance of graphene,” Phys. Rev. B Condens. Matter **87**(11), 115429 (2013). [CrossRef]

**20. **D. Correas-Serrano, J. S. Gomez-Diaz, and A. Alvarez-Melcon, “On the Influence of Spatial Dispersion on the Performance of Graphene-Based Plasmonic Devices,” IEEE Anten. Wirel. Propag. Lett. **13**, 345–348 (2014).

**21. **D. Correas-Serrano, J. S. Gomez-Diaz, J. Perruisseau-Carrier, and A. Alvarez-Melcon, “Spatially Dispersive Graphene Single and Parallel Plate Waveguides: Analysis and Circuit Models,” IEEE Trans. Microw. Theory Tech. **61**(12), 4333–4344 (2013). [CrossRef]

**22. **J. S. Gomez-Diaz, J. Mosig, and A. Alvarez-Melcon, “Effect of Spatial Dispersion on Surface Waves Propagating Along Graphene Sheets,” IEEE Trans. Antenn. Propag. **61**, 3589–3596 (2013).

**23. **A. Fallahi, T. Low, M. Tamagnone, and J. Perruisseau-Carrier, “Nonlocal electromagnetic response of graphene nanostructures,” Phys. Rev. B **91**(12), 121405 (2015). [CrossRef]

**24. **M. Tymchenko, A. Y. Nikitin, and L. Martín-Moreno, “Faraday Rotation Due to Excitation of Magnetoplasmons in Graphene Microribbons,” ACS Nano **7**(11), 9780–9787 (2013). [CrossRef] [PubMed]

**25. **O. Luukkonen, C. Simovski, G. Grant, G. Goussetis, D. Lioubtchenko, A. Raisanen, and S. Tretyakov, “Simple and Accurate Analytical Model of Planar Grids and High-Impedance Surfaces Comprising Metal Strips or Patches,” IEEE Trans. Antenn. Propag. **56**, 1624 (2008).

**26. **E. Forati, G. W. Hanson, A. B. Yakovlev, and A. Alù, “Planar hyperlens based on a modulated graphene monolayer,” Phys. Rev. B **89**(8), 081410 (2014). [CrossRef]

**27. **V. P. Gusynin, S. G. Sharapov, and J. P. Carbotte, “Magneto-optical conductivity in graphene,” J. Phys. Condens. Matter **19**(2), 026222 (2007). [CrossRef]

**28. **I. Crassee, M. Orlita, M. Potemski, A. L. Walter, M. Ostler, T. Seyller, I. Gaponenko, J. Chen, and A. B. Kuzmenko, “Intrinsic Terahertz Plasmons and Magnetoplasmons in Large Scale Monolayer Graphene,” Nano Lett. **12**(5), 2470–2474 (2012). [CrossRef] [PubMed]

**29. **J. S. Gómez-Díaz and J. Perruisseau-Carrier, “Propagation of hybrid transverse magnetic-transverse electric plasmons on magnetically biased graphene sheets,” J. Appl. Phys. **112**(12), 124906 (2012). [CrossRef]

**30. **W. L. Mochan, M. del CastilloMussot, and R. G. Barrera, “Effect of plasma waves on the optical properties of metal-insulator superlattices,” Phys. Rev. B **35**(3), 1088–1098 (1987). [CrossRef] [PubMed]

**31. **M. G. Silveirinha, “Additional boundary condition for the wire medium,” IEEE Trans. Antenn. Propag. **54**, 1766–1780 (2006).

**32. **A. D. Boardman, *Electromagnetic Surface Modes*, John Wiley and Sons, Chichester (1982).

**33. **H. J. Bilow, “Guided Waves on a Planar Tensor Impedance Surface,” IEEE Trans. Antenn. Propag. **51**(10), 2788–2792 (2003). [CrossRef]

**34. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C + + : the art of scientific computing, 3rd edition (Cambridge University Press, 2007).

**35. **L. Chen, C. Zhang, J. Zhou, and L. J. Guo, “Probing the Ultrathin Limit of Hyperbolic Meta-material: Nonlocality Induced Topological Transitions,” in *CLEO**:**2015*, OSA Technical Digest (online) (Optical Society of America, 2015), paper FTu4C.6.

**36. **D. Rodrigo, O. Limaj, D. Janner, D. Etezadi, F. J. García de Abajo, V. Pruneri, and H. Altug, “Mid-infrared plasmonic biosensing with graphene,” Science **349**(6244), 165–168 (2015). [CrossRef] [PubMed]

**37. **A. Alù, “First-principles homogenization theory for periodic metamaterials,” Phys. Rev. B **84**(7), 075153 (2011). [CrossRef]

**38. **J. G. Son, M. Son, K. J. Moon, B. H. Lee, J. M. Myoung, M. S. Strano, M. H. Ham, and C. A. Ross, “Sub-10 nm Graphene Nanoribbon Array Field-Effect Transistors Fabricated by Block Copolymer Lithography,” Adv. Mater. **25**(34), 4723–4728 (2013). [CrossRef] [PubMed]

**39. **L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University Press, 2006).

**40. **N. A. P. Nicorovici, R. C. McPhedran, and L. C. Botten, “Relative local density of states for homogeneous lossy materials,” Phys. B Condens. Matter **405**, 2915 (2010).

**41. **O. Di Stefano, N. Fina, S. Savasta, R. Girlanda, and M. Pieruccini, “Calculation of the local optical density of states in absorbing and gain media,” J. Phys. Condens. Matter **22**(31), 315302 (2010). [CrossRef] [PubMed]