Singular graphene metasurfaces, conductivity gratings realized by periodically suppressing the local doping level of a graphene sheet, were recently proposed to efficiently harvest THz light and couple it to surface plasmons over broad absorption bands, thereby achieving remarkably high field enhancement. However, the large momentum wavevectors thus attained are sensitive to the nonlocal behavior of the underlying electron liquid. Here, we extend the theory of singular graphene metasurfaces to account for the full nonlocal optical response of graphene and discuss the resulting impact on the plasmon resonance spectrum. Finally, we propose a simple local-analogue model that is able to reproduce the effect of nonlocality in local-response calculations by introducing a constant conductivity offset, which could prove a valuable tool in the modeling of more complex experimental graphene-based platforms.
Over the past two decades singular plasmonic structures, such as touching metallic wires and spheres, demonstrated enticing capabilities for controlling light in the subwavelength regime, thanks to their ability to bridge very different length scales, namely, the wavelength of the photon and that of the electron , , . Characterized by features much smaller than their overall size, such as sharp points and regions of adiabatically vanishing thickness, these structures, so far, enabled extreme confinement of electromagnetic fields, with a plethora of far-reaching applications, including the access to quantum regimes of light–matter interactions , , . More recently, extended structures featuring geometrical singularities were investigated in the context of metasurfaces , , which enable larger scattering cross-sections and lower losses, as well as unprecedented tunability and dynamical control of electromagnetic waves , , .
The working principle of singular structures, which was recently shown to be intimately linked to the concept of compactification encountered in high-dimensional field theories , , may be summarized in the following consideration. In a conventional one-dimensional (1D) periodic scattering problem (Figure 1), such as that of a plasmon propagating on a inhomogeneously doped graphene sheet, one can identify two distinct scenarios: hard-boundary scattering, which is often modeled through boundary conditions, commonly results in reflection, and the subsequent quantization of scattered fields into effective Fabry–Pérot modes (Figure 1A); the opposite regime consists of the weak scattering limit, often modeled with WKB-type approaches, whose main effect is the phase change of a largely transmitted wave (Figure 1B), which results in the opening of band-gaps as a result of Bragg scattering from the periodic inhomogeneity. Singular structures constitute a narrow intermediate regime, whereby the scattering process is not abrupt enough to generate significant back-reflection, while not being smooth enough to let the wave be significantly transmitted and interfere. As a result, the wavelength of the excitation becomes increasingly short as it approaches a so-called singular point, which we shall define as a point where the local wavevector effectively diverges. Its group velocity is dramatically reduced, such that the wave never reaches the singularity, and energy is absorbed close to it in the presence of material loss, realizing a remarkable concentration of electromagnetic energy within nanoscale volumes. Recently, graphene-based singular metasurfaces were proposed as a promising platform for the focusing of THz plasmons, as well as for their broadband, tunable plasmonic response to far-field illumination . The plasmonic response of graphene recently demonstrated unprecedented field confinement, concentrating waves which propagate with free-space wavelengths of tens to hundreds of microns down to the atomic scale , , , . In addition, the technological relevance of these THz plasmons for vibrational sensing , , , ,  and high-speed wireless communication , ,  attracted enormous interest in these surface excitations.
However, it was recently shown that the account of nonlocal effects – arising from the quantum nonlocal response of the 2D electron gas – is of paramount importance when the plasmon wavelength becomes comparable to the electronic Fermi wavelength, in order to correctly predict their electromagnetic response , . The nonlocal response of singular metallic structures featuring 3D electron gases was widely studied , primarily via the so-called hydrodynamic model , , which accounts for charge screening at a dielectric–metal interface , , . Alternative theoretical models were also developed in the past, which simplify the account of nonlocal effects in complex plasmonic structures , , , . More recently, nonlocal effects attracted renewed interest and, in particular, due to the sizable impact of quantum mechanical effects in plasmon-enhanced light–matter interactions ,  at the nanoscale, as well as for applications to all-optical signal processing .
In these singular metasurfaces, the nonlocal response of graphene arises from the onset of different types of electronic transitions within the regions of phase space shown in Figure 2. Region 1B constitutes the so-called lossless regime (in the absence of electronic scattering processes). Here, interband transitions are forbidden due to Pauli blocking, and the small plasmon momentum – i.e. k≪kF, where kF is the Fermi wavevector – does not allow for any indirect transitions. Hence, in this regime, the only loss channels for graphene plasmons arise from electronic scattering processes (e.g. with phonons, defects, etc.) , , which are commonly introduced phenomenologically via the so-called relaxation-time approximation . Nevertheless, the incorporation of quantum nonlocal effects is reflected in the reactive (imaginary) component of the conductivity for large plasmon momenta k→→ω/vF. In fact, the divergent character of graphene’s conductivity at the boundary between region 1B and 1A constitutes a main detrimental effect for the realization of conductivity singularities in graphene. Region 1A suffers from the onset of Landau damping, which arises due to the matching between the phase velocities of the electrons and of the plasmons. This has the effect of dramatically enhancing the loss. Similarly, region 2A is affected by additional intraband channels, which become accessible once the plasmon momentum k>kF. Finally, indirect (region 2B) and direct (region 3B) interband transitions occur once the plasmon energy ħω>2EF−ħvFk and ħω>2EF+ħvFk, respectively.
Because of the extreme values of plasmon momenta to which a singular structure can couple incident photons to, a rigorous account of the momentum dependence of the optical response of these metasurfaces is pivotal. In this work, we explore the nonlocal behavior of plasmons in singular graphene metasurfaces and show that these systems are able to probe the strong nonlocal response of 2D electron gases by coupling far-field radiation to deeply subwavelength plasmon modes. By means of a nonlocal mode-matching technique , supported by numerical calculations, as well as a phenomenological local-analogue model, we unravel the physics underpinning the onset of nonlocality in these metasurfaces. We believe that our method constitutes a valuable tool for incorporating nonlocal effects in complex metasurface setups and may be employed as an alternative approach to fully nonlocal conductivity models.
Nonlocal effects in plasmonics manifest themselves when the plasmon wavelength approaches the typical electronic wavelength λF in a material. In this regime, the spatial variation of the electric field E(x) is sufficiently abrupt to sample the underlying inhomogeneity of the electron gas, so that the constitutive relation for the surface current density can be written as
and thus can no longer be approximated assuming a spatial dependence of the conductivity of the form σ(x−x′, ω)= σ(ω)δ(x−x′), where δ(x) is the Dirac delta function.
However, when the structuring of a THz metasurface is performed over scales much larger than the Fermi’s wavelength (L≫λF), a separation of length scales can be assumed. Hence, we can write, under the adiabatic approximation:
where ζ(x′) is a dimensionless variable, which describes the spatial modulation of the conductivity of graphene , , the latter depending monotonically on the local doping level of graphene. This has the desirable property of being actively tunable (e.g. electrostatically, chemically, or optically). In this work, we assume that a periodic conductivity modulation is applied, which, for simplicity and definiteness, is herein assumed to be of the form ζ(x)=1+ζ1cos(gx), where L=2π/g is the period of the 1D metasurface and g the reciprocal lattice vector associated with the same.
Using Bloch’s theorem and expanding the Bloch modes of the in-plane electric field and the surface current as a Fourier series, one may write
and a simple relation between the Fourier amplitudes of the electric field and the surface current hereby takes the form
which is accurate as long as the reciprocal lattice vector of the metasurface satisfies g≪kF. For concreteness, the nonlocal conductivity model  is described in Appendix 1.
The main effect of nonlocality in graphene is to oppose the formation of a singularity by increasing the conductivity probed by large-momentum Fourier components. In Figure 3, we plot the transmission spectra under plane wave illumination at normal incidence (k=0) for different modulation strengths Δ=−log10(1−ζ1), corresponding to the number of orders of magnitude by which the conductivity is suppressed at a singular point. We assume an average Fermi level EF=0.4 eV, a conductivity grating period L=5 μm, and a mobility μm=104 cm2/(V s) resulting in an electron scattering time where the Fermi velocity vF=9.5×105 m/s  is assumed. Our results are obtained via the nonlocal mode-matching method outlined above; these were benchmarked, in the local-response limit, against finite-element method (FEM) numerical calculations using a commercially available package (COMSOL Multiphysics). For weak conductivity modulation, i.e. far from the singular limit (Figure 3A), the local and nonlocal spectra are effectively equivalent. In this limit, only momentum states well below the Landau damping regime k≈ω/vF are populated, so that the metasurface can be accurately described via a local Drude-type conductivity model where γ=τ−1. As we increase the modulation strength to 99.9% of the average value (Figure 3B, Δ=2), the local and nonlocal spectra start deviating, the latter exhibiting a clear blueshift, which is a consequence of nonlocality (see, e.g. Ref. ), as plasmon resonance frequencies ω∝σ [see dispersion relation, Eq. (5)], and nonlocal effects lead to an increase in conductivity. Finally, for Δ=3 (Figure 3), nonlocality becomes a dominant effect, which effectively saturates the plasmonic spectrum, opposing any further merging of the plasmon resonances.
For completeness we show in Figure 4 the plasmonic band structure of our metasurfaces over the rest of the Brillouin zone, where a few additional effects are present. In order to visualize the bands, we plot in log-scale the absolute value of the reflection coefficient, which was color saturated in order to allow both propagating and evanescent modes to be identifiable. In the non-singular regime (Figure 4A), plasmonic band gaps resulting from the periodic modulation are clearly visible at k=π/L, whereas the bands are degenerate at k=0 due to the inversion symmetry of the modulation. By contrast, as the modulation strength is increased (Figure 4B), the bands flatten as a result of the stronger Bragg scattering, so that in the singular limit, they become effectively indistinguishable. In this regime, plasmons are dramatically slowed down. However as in the previous case, the introduction of nonlocality saturates the merging of the plasmon spectrum, opposing the flattening of the bands. However, in this case, we clearly see how nonlocality has the additional effect of broadening the reflection spectra dramatically, as a result of the losses introduced by Landau damping.
The account of nonlocality can be somewhat demanding in the modeling of more complex experimental setups. Consequently, local-analogue models, which are able to incorporate the effects of nonlocality in a local simulation are valuable tools for the theoretical modeling of plasmonic systems. Here, we propose a simple local-analogue model, which can accurately reproduce the results of the fully nonlocal calculation carried out above. Local-analogue models were originally proposed for metallic plasmonic systems  in order to capture nonlocal effects under the framework of the hydrodynamic model of the free-electron gas at the interface between nearly touching metallic structures. In that context, the effect of nonlocality is the inward shift of the induced charges, i.e. away from the metallic surface and into the bulk, thereby effectively widening the gap between the components of the dimer (e.g. metallic cylinders or spheres). Consequently, the substitution of a thin metallic layer by an effective dielectric one was able to accurately reproduce the optical response of such nearly touching metallic structures.
Conversely, the type of singular structure described in this work entails the inverse effect: since the conductivity is strongly enhanced as k→ω/vF, the effect of nonlocality is to smear out the singularity by effectively saturating the local conductivity to a minimum level σs dictated, qualitatively, by the condition k(σs)≈ω/vF, i.e. when the plasmon wavelength λp→λF, and Landau damping opposes any further confinement of the plasmonic field. The quasi-static dispersion relation of graphene plasmons is read as :
where σ≡σ(k, ω) and σ≡σ(ω)=σ(k→0, ω) in the nonlocal and local cases, respectively. Herein, we set ε1,2=1 (for simplicity alone). Moreover, we can then substitute the wavevector k=βω/vF, where β is a phenomenological factor of order ~1, which quantifies the fraction of electron momentum to which the plasmon can couple before saturating (which is exactly one if momentum saturation occurs exactly at the electron momentum). In this fashion, we thus obtain the saturation value for the conductivity, σs=2iε0vF/β. In Figure 5, we add a positive surface conductivity offset
in a local FEM calculation, where the factor in the first square bracket is responsible for the smearing of the imaginary part of the surface conductivity, whereas the second ensures that the loss tangent is preserved upon the conductivity offset. It is worth remarking that, as our model hinges on the relation between the plasmonic and electronic momentum, the linear dependence of the latter on frequency implies that the conductivity offset is frequency dependent.
For β=1, the agreement between the previous nonlocal result (Figure 3C) and the spectrum obtained using the local-analogue model is only qualitative. However, as the figure plainly shows, by choosing β≃1.29, this simple model is able to reproduce the entire transmission spectrum with remarkable accuracy, hereby validating the physical assumptions behind our local-analogue model. While insightful, our semiclassical theory does not provide a quantitative evaluation of the saturation parameter β. However, it can shine further insight into our nonlocal description: in fact, given that β>1, meaning that the saturation momentum surpasses the electron momentum, this result is worthy of a closer inspection. In the most singular regime, the assumption leading to the first-order Fourier expansion in Eq. (2) may lose accuracy for larger wavevectors, resulting in an underestimation of the extent of nonlocal effects in this extreme regime. In this sense, comparisons with future, fully quantum-mechanical investigations would prove extremely useful in providing an exact evaluation of the saturation momentum. Nevertheless, our local-analogue model offers a useful and intuitive method for the incorporation of nonlocal effects in the future modeling of complex metasurfaces based on 2D materials.
In this work, we presented a theoretical description of nonlocal effects in singular graphene metasurfaces. By calculating the transmission spectra under plane wave illumination, as well as the plasmon band structure, we demonstrated how such conductivity gratings are able to probe the nonlocal response of graphene. Furthermore, we discussed the consequent limitations imposed by nonlocality to the field confinement and spectral degeneracy induced by the singularity, which is effectively smeared out by the increased conductivity probed by large plasmon wavevectors. Finally, we proposed a simple local-analogue model, which is able to reproduce the effects of nonlocality by means of an effective surface conductivity offset, which saturates the plasmon wavevector to the electronic one. An analogous effective model could also be devised to account for nonlocality in metallic surfaces with singular points, extending previous work on local-analogue models for metal nanostructures . This would find particular interest in the context of the recent advances in the fabrication of ultrathin metals , .
To conclude, singular graphene metasurfaces constitute a platform for probing nonlocality in graphene with far-field measurements. Our results form the basis for a quantitative account of nonlocality in these 2D systems, and should be valuable for guiding future experimental efforts, as well as incorporating fully quantum-mechanical theoretical investigations into effective local descriptions.
We would like to thank A.I. Fernández-Domínguez and F. Koppens for useful discussions. E.G. was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the EPSRC (EP/L015579/1). P.A.H. acknowledges funding from Fundação para a Ciência e a Tecnologia and Instituto de Telecomunicações under projects CEECIND/03866/2017 and UID/EEA/50008/2019. J.B.P. acknowledges funding from the Gordon and Betty Moore Foundation. N.A.M. is a VILLUM Investigator supported by VILLUM FONDEN (grant No. 16498). The Center for Nano Optics is financially supported by the University of Southern Denmark (SDU 2020 funding). The Center for Nanostructured Graphene is sponsored by the Danish National Research Foundation (Project No. DNRF103).
Nonlocal conductivity model
The nonlocal conductivity of graphene can be written in terms of graphene’s 2D polarizability as 
where Pγ(k, ω) is the 2D density-density response function (or 2D polarizability) in the relaxation-time approximation (which incorporates a finite plasmon lifetime while preserving electron number density , ). The 2D polarizability in the relaxation-time approximation is given by , 
where P(k, ω) denotes the zero-temperature density-density response function in the four regions outlined in Figure 2, which may be written as:
where the constant and the auxiliary functions:
 Luo Y, Pendry JB, Aubry A. Surface plasmons and singularities. Nano Lett 2010;10:4186–91. Search in Google Scholar
 Aubry A, Lei DY, Fernández-Domnguez AI, Sonnefraud Y, Maier SA, Pendry JB. Plasmonic light-harvesting devices over the whole visible spectrum. Nano Lett 2010;10:2574–9. Search in Google Scholar
 Estakhri NM, Alù A. Physics of unbounded, broadband absorption/gain efficiency in plasmonic nanoparticles. Phys Rev B 2013;87:205418. Search in Google Scholar
 Chikkaraddy R, De Nijs B, Benz F, Barrow SJ, Scherman OA, Rosta E. Single-molecule strong coupling at room temperature in plasmonic nanocavities. Nature 2016;535:127. Search in Google Scholar
 Zhu W, Esteban R, Borisov AG, Baumberg JJ, Nordlander P, Lezec HJ. Quantum mechanical effects in plasmonic structures with subnanometre gaps. Nat Commun 2016;7:11495. Search in Google Scholar
 Fernández-Domínguez AI, Bozhevolnyi SI, Mortensen NA. Plasmon-enhanced generation of nonclassical light. ACS Photonics 2018;5:3447–51. Search in Google Scholar
 Galiffi E, Pendry JB, Huidobro PA. Broadband tunable THz absorption with singular graphene metasurfaces. ACS Nano 2018;12:1006–13. Search in Google Scholar
 Yang F, Huidobro PA, Pendry JB. Transformation optics approach to singular metasurfaces. Phys Rev B 2018;98:125409. Search in Google Scholar
 Kildishev AV, Boltasseva A, Shalaev VM. Planar photonics with metasurfaces. Science 2013;339:1232009. Search in Google Scholar
 Yu N, Capasso F. Flat optics with designer metasurfaces. Nat Mater 2014;13:139. Search in Google Scholar
 Shaltout AM, Shalaev VM, Brongersma ML. Spatiotemporal light control with active metasurfaces. Science 2019;364:eaat3100. Search in Google Scholar
 Pendry JB, Huidobro PA, Luo Y, Galiffi E. Compacted dimensions and singular plasmonic surfaces. Science 2017;358:915–7. Search in Google Scholar
 Galiffi E, Pendry JB, Huidobro PA. Singular graphene metasurfaces. EPJ Appl Metamater 2019;6:10. Search in Google Scholar
 Iranzo DA, Nanot S, Dias EJC, Epstein I, Peng C, Efetov DK. Probing the ultimate plasmon confinement limits with a van der Waals heterostructure. Science 2018;360:291–5. Search in Google Scholar
 Lundeberg MB, Gao Y, Asgari R, Tan C, Van Duppen B, Autore M. Tuning quantum nonlocal effects in graphene plasmonics. Science 2017;357:187–91. Search in Google Scholar
 Dias EJC, Iranzo DA, Gonçalves PAD, Hajati Y, Bludov YV, Jauho A-P. Probing nonlocal effects in metals with graphene plasmons. Phys Rev B 2018;97:245405. Search in Google Scholar
 Ni GX, McLeod AS, Sun Z, Wang L, Xiong L, Post KW. Fundamental limits to graphene plasmonics. Nature 2018;557:530–3. Search in Google Scholar
 Rodrigo D, Limaj O, Janner D, et al. Mid-infrared plasmonic biosensing with graphene. Science 2015;349:165–8. Search in Google Scholar
 Li Y, Yan H, Farmer DB, Meng X, Zhu W, Osgood RM. Graphene plasmon enhanced vibrational sensing of surface-adsorbed layers. Nano Lett 2014;14:1573–7. Search in Google Scholar
 Gonçalves PAD, Peres NMR. An introduction to graphene plasmonics, 1st ed. Singapore: World Scientific, 2016. Search in Google Scholar
 Grigorenko AN, Polini M, Novoselov KS. Graphene plasmonics. Nat Photonics 2012;6:749. Search in Google Scholar
 Low T, Avouris P. Graphene plasmonics for terahertz to mid-infrared applications. ACS Nano 2014;8:1086–101. Search in Google Scholar
 Vicarelli L, Vitiello MS, Coquillat D, Lombardo A, Ferrari AC, Knap W. Graphene field-effect transistors as room-temperature terahertz detectors. Nat Mater 2012;11:865–71. Search in Google Scholar
 Poumirol JM, Liu PQ, Slipchenko TM, Nikitin AY, Martín-Moreno L, Faist J. Electrically controlled terahertz magneto-optical phenomena in continuous and patterned graphene. Nat Commun 2017;8:14626. Search in Google Scholar
 Yang X, Vorobiev A, Generalov A, Andersson MA, Stake J. A flexible graphene terahertz detector. Appl Phys Lett 2017;111:021102. Search in Google Scholar
 Raza S, Bozhevolnyi SI, Wubs M, Mortensen NA. Nonlocal optical response in metallic nanostructures. J Phys Cond Mater 2015;27:183204. Search in Google Scholar
 Boardman AD. Electromagnetic surface modes. Chichester, New York, John Wiley & Sons, 1982. Search in Google Scholar
 Schwartz C, Schaich WL. Hydrodynamic models of surface plasmons. Phys Rev B 1982;26:7008–11. Search in Google Scholar
 García de Abajo FJ. Nonlocal effects in the plasmons of strongly interacting nanoparticles, dimers, and waveguides. J Phys Chem C 2008;112:17983–7. Search in Google Scholar
 Toscano G, Raza S, Yan W, Jeppesen C, Xiao S, Wubs M. Nonlocal response in plasmonic waveguiding with extreme light confinement. Nanophotonics 2013;2:161–6. Search in Google Scholar
 Yang F, Wang YT, Huidobro PA, Pendry JB. Nonlocal effects in singular plasmonic metasurfaces. Phys Rev B 2019;99:165423. Search in Google Scholar
 Luo Y, Fernández-Domínguez AI, Wiener A, Maier SA, Pendry JB. Surface plasmons and nonlocality: a simple model. Phys Rev Lett 2013;111:093901. Search in Google Scholar
 Mortensen NA, Raza S, Wubs M, Søndergaard T, Bozhevolnyi SI. A generalized nonlocal optical response theory for plasmonic nanostructures. Nat Commun 2014;5:3809. Search in Google Scholar
 Yan W, Wubs M, Mortensen NA. Projected dipole model for quantum plasmonics. Phys Rev Lett 2015;115:137403. Search in Google Scholar
 Christensen T, Yan W, Jauho AP, Soljačić M, Mortensen NA. Quantum corrections in nanoplasmonics: Shape, scale, and material. Phys Rev Lett 2017;118:157402. Search in Google Scholar
 Gonçalves PAD, Christensen T, Rivera N, Jauho AP, Mortensen NA, Soljačić M. Plasmon-emitter interactions at the nanoscale. 2019;arXiv:1904.09279. Search in Google Scholar
 Kwon H, Sounas D, Cordaro A, Polman A, Alù A. Nonlocal metasurfaces for optical signal processing. Phys Rev Lett 2018;121:173004. Search in Google Scholar
 Yan H, Low T, Zhu W, Wu Y, Freitag M, Li X. Damping pathways of mid-infrared plasmons in graphene nanostructures. Nat Photonics 2013;7:394–9. Search in Google Scholar
 Slipchenko TM, Nesterov ML, Martín-Moreno L, Nikitin AY. Analytical solution for the diffraction of an electromagnetic wave by a graphene grating. J Opt 2013;15:114008. Search in Google Scholar
 Huidobro PA, Kraft M, Kun R, Maier SA, Pendry JB. Graphene, plasmons and transformation optics. J Opt 2016;18:044024. Search in Google Scholar
 Fan Y, Shen NH, Koschny T, Soukoulis CM. Tunable terahertz meta-surface with graphene cut-wires. ACS Photonics 2015;2:151–6. Search in Google Scholar
 El-Fattah Abd ZM, Mkhitaryan V, Brede J, et al. Plasmonics in atomically thin crystalline silver films. ACS Nano 2019;13:7771–9. Search in Google Scholar
 Maniyara RA, Rodrigo D, Yu R, Canet-Ferrer J, Ghosh DS, Yongsunthon R. Tunable plasmons in ultrathin metal films. Nat Photonics 2019;13:328. Search in Google Scholar
 Mermin ND. Lindhard dielectric function in the relaxation-time approximation. Phys Rev B Mar 1970;1:2362–3. Search in Google Scholar
©2020 Emanuele Galiffi et al., published by De Gruyter, Berlin/Boston
This work is licensed under the Creative Commons Attribution 4.0 Public License.