1 Introduction

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 [1], [2], [3]. 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 [4], [5], [6]. More recently, extended structures featuring geometrical singularities were investigated in the context of metasurfaces [7], [8], which enable larger scattering cross-sections and lower losses, as well as unprecedented tunability and dynamical control of electromagnetic waves [9], [10], [11].

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 [12], [13], 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 [7]. 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 [14], [15], [16], [17]. In addition, the technological relevance of these THz plasmons for vibrational sensing [18], [19], [20], [21], [22] and high-speed wireless communication [23], [24], [25] attracted enormous interest in these surface excitations.

Figure 1:

Traveling towards a geometrical singularity, a wave is unable to be reflected or transmitted, becoming increasingly compressed.

(A) The in-plane scattering of an electromagnetic wave in a periodic system, e.g. a plasmon propagating along a periodically modulated conductive surface is typically dominated by reflection at hard boundaries or transmission through soft boundaries, leading to discrete Fabry–Pérot modes or Bloch waves, respectively. (B) At a singular boundary, both transmission and reflection channels are virtually inaccessible, and the only available path for a wave is to shrink its wavelength and concentrate its energy as it travels toward the singular point.

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 [15], [16]. The nonlocal response of singular metallic structures featuring 3D electron gases was widely studied [26], primarily via the so-called hydrodynamic model [27], [28], which accounts for charge screening at a dielectric–metal interface [29], [30], [31]. Alternative theoretical models were also developed in the past, which simplify the account of nonlocal effects in complex plasmonic structures [32], [33], [34], [35]. 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 [16], [36] at the nanoscale, as well as for applications to all-optical signal processing [37].

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≪k_F, where k_F 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.) [17], [38], which are commonly introduced phenomenologically via the so-called relaxation-time approximation [20]. Nevertheless, the incorporation of quantum nonlocal effects is reflected in the reactive (imaginary) component of the conductivity for large plasmon momenta k→→ω/v_F. 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>k_F. Finally, indirect (region 2B) and direct (region 3B) interband transitions occur once the plasmon energy ħω>2E_F−ħv_Fk and ħω>2E_F+ħv_Fk, respectively.

Figure 2:

Electronic contributions to the graphene conductivity in different regions of phase space [20].

Region 1B (k<ω/v_F, k<2k_F−ω/v_F) of phase space is protected from Landau damping arising from both interband and intraband transitions. The lossy (shaded) regions are 1A (ω/v_F<k<2k_F−ω/v_F) and 2A (ω/v_F<k<2k_F+ω/v_F, k>2k_F−ω), dominated by Landau damping resulting from intraband transitions, and 2B (2k_F−ω/v_F<k<ω/v_F, ω/v_F<2k_F+k) and 3B (k<ω/v_F−2k_F) dominated by indirect and direct interband transitions, 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 [20], 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.

2 Methods

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 [39], [40], 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+ζ₁cos(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≪k_F. For concreteness, the nonlocal conductivity model [20] is described in Appendix 1.

3 Results

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 Δ=−log₁₀(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 E_F=0.4 eV, a conductivity grating period L=5 μm, and a mobility μ_m=10⁴ cm²/(V s) resulting in an electron scattering time τ=μmEF/(vF2e)≈0.44  ps, where the Fermi velocity v_F=9.5×10⁵ m/s [41] 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≈ω/v_F are populated, so that the metasurface can be accurately described via a local Drude-type conductivity model σD(ω)=e2πℏ2EF(γ−iω), where γ=τ⁻¹. 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. [20]), 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.

Figure 3:

Nonlocality leads to saturation of the density of states in a singular graphene metasurface.

Local (red) and nonlocal (blue) transmittance spectra for plane wave illumination through the graphene metasurface at normal incidence, obtained with the mode-matching (continuous lines) and finite-element method (dots) for three increasingly singular metasurfaces corresponding to Δ=1 (A), Δ=2 (B), and Δ=3 (C), respectively. The nonlocal contribution, which is negligible away from the singular regime, becomes dominant as the singular limit is approached, opposing the merging of surface plasmon modes.

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.

Figure 4:

(Left-to-right) Plasmonic band structures for non-singular, singular (local) and singular (nonlocal) graphene metasurfaces.

Plasmonic band structure of the non-singular (Δ=1) (A) and singular (Δ=3) (B, C) metasurfaces, visualized by plotting the logarithm of the absolute value of the reflection coefficient. Note that to the right of the light line, the incoming waves are evanescent. Local (A, B) and nonlocal (C) spectra differ significantly for the singular case Δ=3 only. (A) In the non-singular regime, band gaps are clearly visible at the edge of the Brillouin zone, with the respective upper bands showing significantly larger broadening compared to the lower ones. (B) In the singular limit, local calculations predict that the band above of each band gap becomes indistinguishable from the lower band of the gap above, effectively realizing a series of flat bands, with extremely low group velocities. However, the onset of Landau damping in the nonlocal case (C) greatly broadens these bands, in addition to saturating their compression as in Figure 3. In addition, note that the bands are effectively degenerate at k=0, due to the inversion symmetry of the modulation.

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 [32] 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→ω/v_F, 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)≈ω/v_F, 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 [20]:

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=βω/v_F, 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ε₀v_F/β. In Figure 5, we add a positive surface conductivity offset

Figure 5:

Local (red), nonlocal (blue line), and local-analogue (green triangles for β=1.29 and gray dashed line for β=1) transmittance spectra of the singular (Δ=3) graphene conductivity grating.

The inset shows how a local-analogue metasurface can be obtained by saturating the conductivity of graphene near the value, which causes the local plasmon dispersion to cross the electron dispersion ω=v_Fk, a regime dominated by Landau damping.

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.

4 Conclusions

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 [32]. This would find particular interest in the context of the recent advances in the fabrication of ultrathin metals [42], [43].

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.