Home | About | Journals | Submit | Contact Us | Français |

**|**Scientific Reports**|**PMC5396075

Formats

Article sections

Authors

Related links

Sci Rep. 2015; 5: 9535.

Published online 2015 April 9. doi: 10.1038/srep09535

PMCID: PMC5396075

Weihua Wang,^{1,}^{2} Thomas Christensen,^{1,}^{2} Antti-Pekka Jauho,^{1,}^{3} Kristian S. Thygesen,^{1,}^{4} Martijn Wubs,^{1,}^{2} and N. Asger Mortensen^{a,}^{1,}^{2}

Received 2014 December 2; Accepted 2015 March 5.

Copyright © 2015, Macmillan Publishers Limited. All rights reserved

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder in order to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

In classical electrodynamics, nanostructured graphene is commonly modeled by the computationally demanding problem of a three-dimensional conducting film of atomic-scale thickness. Here, we propose an efficient alternative two-dimensional electrostatic approach where all calculation procedures are restricted to the graphene sheet. Furthermore, to explore possible quantum effects, we perform tight-binding calculations, adopting a random-phase approximation. We investigate multiple plasmon modes in 20 nm equilateral triangles of graphene, treating the optical response classically as well as quantum mechanically. Compared to the classical plasmonic spectrum which is “blind” to the edge termination, we find that the quantum plasmon frequencies exhibit blueshifts in the case of armchair edge termination of the underlying atomic lattice, while redshifts are found for zigzag edges. Furthermore, we find spectral features in the zigzag case which are associated with electronic edge states not present for armchair termination. Merging pairs of triangles into dimers, plasmon hybridization leads to energy splitting that appears strongest in classical calculations while splitting is lower for armchair edges and even more reduced for zigzag edges. Our various results illustrate a surprising phenomenon: Even 20 nm large graphene structures clearly exhibit quantum plasmonic features due to atomic-scale details in the edge termination.

The collective excitations of conduction electrons in noble metals have been of great interest for a very long time. These excitations known as plasmons play an important role in the optical properties of metals. Through strong plasmon-photon interactions, metals can support important phenomena, such as focusing beyond the diffraction limit^{1}, squeezing the light down to nanoscale^{2}, and large local field enhancement^{3}. Due to these features, plasmons in metals give rise to various potential applications, and especially form a bridge between the worlds of photonics and electronics which commonly work at different length scales^{4}. Developments in nanofabrication technology have stimulated a series of plasmon-based devices like waveguides^{5}, filters^{6}, switches^{7}, and modulators^{8}. In many respects, plasmonic devices open a door to a better performance in speed and size, holding potential for faster dynamics than electronic devices while still having a smaller size footprint than the common all-dielectric photonic devices. However, the inherent Joule loss in metals severely hampers many practical applications of plasmonics^{9}. Alternatively, attempts have already been made to study plasmonics in materials other than metals^{10}, for example doped semiconductors^{11} and superconductors^{12}^{,13}.

Graphene and other low-dimensional crystals are now emerging as interesting materials for exciting science and technology^{14}. Here we study the plasmonic properties of graphene flakes. In its pristine form graphene is a semimetal, but with appropriate doping it is emerging as a promising plasmonic material as well^{15}^{,16}^{,17}^{,18}^{,19}. The graphene plasmons are non-radiating, but with a momentum mismatch to free-space radiation that can be overcome with the aid of e.g. grating approaches^{20}. The charge carriers in graphene obey linear energy dispersion at lower energies close to the Dirac points, thus resembling the linear dispersion of photons^{21}^{,22}^{,23}. Experimental investigations of carrier transport show that the mobility limited by impurity scattering can exceed 15.000 cm^{2}/Vs at room temperature^{21}, which gives the intrinsic loss in graphene one order of magnitude less than the noble metals. Despite relaxation due to phonon scattering^{24}^{,25}, graphene achieves superior plasmonic performance in propagation length and field enhancement^{26}^{,27}. The carrier density in graphene may be adjusted by electrostatic gating, which results in actively tunable plasmons beyond structural variations in metals, as has already been demonstrated experimentally^{28}^{,29}. With the typical doping levels, the plasmonic response is generally in the terahertz (THz) to mid-infrared frequency range, thus allowing new progress in THz technology^{30}. As an example, graphene waveguides (with sub-wavelength width) pave a promising way to realize ultra-compact THz devices where bends and splitters do not bring any significant loss^{31}.

Because of these attractive plasmonic properties, it is worth to comprehensively study the optical properties of graphene. Here the fundamental quantity is the dielectric function. For graphene systems, the dielectric function can be obtained within the framework of linear-response theory and the random-phase approximation (RPA)^{15}^{,16}^{,32}. For infinite graphene sheets, the derived two-dimensional (2D) dielectric function *ε*(*q*,*ω*) is a function of both frequency *ω* and momentum *q*. This is different from common three-dimensional (3D) photonic materials which are usually well-described by frequency-dependent functions, while spatial dispersion is negligible for good dielectrics and most metals (beyond the nanoscale). Two common approximations in the modelling of graphene structures are to adopt the local-response approximation (applying the small-*q* limit) and to model graphene as a very thin conducting film, yet preserving its 3D representation^{33}^{,34}. Using dielectric functions so obtained, one can solve Maxwell's equations for arbitrarily shaped flakes of nanostructured graphene. For very small flakes of characteristic dimension *R* (*R* ~ *λ*_{F} with *λ*_{F} ~ 10 nm being the Fermi wavelength corresponding to a Fermi level of = 0.4 eV), the common assumption is jeopardized and nonlocal response turns important for far-field optical properties. Obviously, near-field properties may be influenced too, e.g. for dimers of sufficiently large structures, where holds, while instead the tiny gap () promotes nonlocal effects^{35}. For optical excitation in the near field^{17}^{,36}, nonlocal response can also be important for large plasmonic structures provided that the distance to the emitter is comparable to the nonlocal characteristic length scale^{37}. In this regime, both semiclassical hydrodynamic^{38}^{,39} and full quantum approaches have been proposed^{40}^{,41}, similar to those recently developed for metals^{35}^{,42}. While previous studies have mainly focused on the optically bright dipole mode, here we will illustrate that structured graphene is also rich on higher-order modes. Although the latter are typically not excited by far-field radiation, they may be probed by near-field optical spectroscopy and/or electron energy loss spectroscopy (EELS).

In this article, we study plasmon properties in individual graphene nanostructures and in dimers of such structures by means of both classical and quantum methods. In particular, we consider triangles of graphene and bow-tie structures formed by such triangles, while our methods can also be applied to other geometries (as we show in the Supplementary Information). Here, we focus on the plasmonic aspects due to doping with a significant number of electrons, while there are also appealing aspects in the single-electron doping regime^{43}. We will emphasize dimers formed by electrically mutually disconnected graphene islands, while graphene dimers connected by quantum junctions^{44} and extended complimentary structures (electronic^{45} and plasmonic^{25}^{,46} periodic anti-dot arrays) also represent interesting geometries and regimes. As a key element, we consider the eigenmodal properties of the plasmonic excitations, thus extending the current state of analysis beyond dipolar excitations, as relevant for near-field plasmonic interaction with e.g. emitters or fast electrons.

In our classical electrodynamical considerations, we treat the nanostructures as 2D materials characterized by a smooth surface conductivity (employing the sheet conductivity derived for bulk graphene), and formulate a closed-form eigenvalue problem on a 2D domain. Numerical solutions in arbitrarily shaped geometries are enabled by finite-element calculations. By its nature, this classical approach neglects the atomic details of the graphene flake. Some aspects e.g. of zigzag termination can be effectively accounted for by additional conductive channels^{39}, though we will not pursue this scheme for our classical calculations presently.

In our quantum treatment, we employ a tight-binding description^{40}^{,41} to account for the actual position of all atoms in the flake and in particular the edge atoms which have the possibility for either armchair or zigzag configurations (Other edge structures can arise from the mixture of these two configurations, but they will not be discussed here). In both the classical and the quantum calculations, multiple plasmon modes are extracted including dipole, multipole, and breathing modes. Their hybridized counterparts in bow-tie nanostructures are also discussed. We show that plasmon excitations and hybridizations are extremely sensitive to the electronic edge effects. This illustrates how quantum plasmonics can manifest itself in graphene structures with dimensions much exceeding the length scales for nonlocal response in individual noble-metal nanoparticles^{35}.

Modern computational electromagnetics is commonly optimized to explore the interaction of radiation with matter in a three-dimensional space, so that two-dimensional material problems are typically not efficiently addressed with existing numerical schemes. For example, a pragmatic approach is to simply mimic the atomically thin graphene layer with a homogenous dielectric film of a finite, yet small thickness *t*_{g}. This assumed 3D film has an effective bulk permittivity, *ε*(*ω*) = *ε*_{0} + *iσ*(*ω*)/(*ωt*_{g}), where *ε*_{0} is the vacuum permability and *σ*(*ω*) denotes the surface conductivity as obtained from e.g. the local-response limit of the RPA^{33}^{,34}. Evidently, the artificial thickness *t*_{g} should be chosen sufficiently small compared with all other characteristic and physical dimensions, yet sufficiently large that meshing hopefully stays computationally feasible and the numerical problem remains tractable. Optimizing this thickness tradeoff does not necessarily give an efficient method. An even more critical issue is that there are no formal proofs, at least to the best of our knowledge, that numerically computed fields (in particular near fields relevant for LDOS or emitter dynamics near a graphene structure^{17}^{,36}) would necessarily converge to physically meaningful quantities in the limit *t*_{g} → 0. Alternatively, in nanostructures with high symmetry, e.g. in ribbons^{47}^{,48} or disks^{39}^{,49}, one may take advantage of modal expansion methods^{39}^{,47}^{,48}^{,49} – which, however, is not an appealing choice for more general structures, where limited analytical progress is possible. In the following, we develop a 2D finite-element approach to efficiently solve the electromagnetic problem self-consistently for graphene in terms of the electric potential and induced charge in general structural configurations.

With the typical sub-eV doping levels, plasmonic resonances typically occur in the mid-infrared regime. The associated free-space wavelength (~10 *μ*m) is then much larger than the geometrical extent of the hosting graphene nanostructures (~10–100 nm). For such problems the electrostatic approximation is excellent. As a computationally very attractive consequence, the electric and constitutive response are governed by two coupled scalar equations for the potential and the induced density *ρ*. In particular, we note that the total potential () is governed by Coulomb's law

where ^{ext}() denotes the external potential, *L* is an auxiliary quantity such as the feature length of the structure which makes the surface integral dimensionless, *ρ*(′) the induced surface charge density, and *ε*_{s} = (*ε*_{above} + *ε*_{below})/2 the averaged dielectric constant of the medium above and below graphene. For simplicity, we only consider freely suspended graphene, so we will use *ε*_{s} = *ε*_{0} throughout the remaining part of the paper. The other scalar equation is obtained by inserting the constitutive equation *J*_{2D} = −*σ*(*ω*)_{2D}() into the continuity equation *iωρ*() = _{2D}·*J*_{2D}(), which for ** r** restricted to the plane of the graphene structure gives

with the 2D Laplace operator. Equation (2) is solved subject to the assumption of charge neutrality, i.e. , implying that on the boundary of the domain, with denoting the in-plane surface normal. The density *ρ* in (2) is restricted to the graphene plane. It may be obtained from a closed-form equation by eliminating the potential in (2) with the help of (1) (see Methods for additional details)^{50}. Once *ρ* within the graphene plane is thus obtained, the potential in the entire space can be evaluated via (1).

Within the framework of the finite-element method (FEM), Eqs. (1) and (2) can both be recast as matrix equations. Concretely, by denoting the FEM-discretized potentials and induced charge densities by vectors, we find the equations ** = **^{ext} + (4*πε _{s}L*)

where **A** and **B** are geometry-dependent square matrices, respectively representing the Coulomb integral in Eq. (1) and the Laplacian in Eq. (2) (see the Methods section below for additional details), while *f*(*ω*) = *iσ*(*ω*)/(4*πε _{s}Lω*) is a geometry-independent scalar

In a quantum mechanical formalism, there are two key computational components: (i) electronic band structure, and (ii) determination of response functions. The graphene *π* and *π** bands (valence and conduction bands respectively) originating from the carbon orbitals are well separated in energy from the four *σ* bands arising from *sp*^{2} hybridization. The dynamics of low-energy excitations in graphene is well-described by inclusion of just the *π* bands, which can be determined by a simple tight-binding model in a nearest-neighbor approximation^{51}^{,52}. Specifically, a graphene nanostructure with *N* carbon atoms results in an *N* × *N* matrix representation of the tight-binding Hamiltonian with elements determined by the orbital hopping integral. A direct diagonalization of the Hamiltonian yields *N* eigenvalues and eigenvectors, corresponding to the electronic energy levels and the wave functions, respectively. The non-interacting density response function, or polarizability matrix , is then built from the electronic states whose elements are given by^{15}^{,16}^{,32}

where denotes the Fermi–Dirac distribution function associated with the state with energy and wave function (*l* labels each of the carbon atoms), while *k*_{B} and are Boltzmann's and Planck's constants, respectively. The factor 2 accounts for spin degeneracy in the absence of a static magnetic field with no Zeeman splitting. In both classical (also called semi-classical due to the conductivity including Fermi–Dirac distribution function) and quantum calculations, states are populated in accordance with a Fermi level of and a temperature *T* = 300 K corresponding to a thermal energy of . We phenomenologically account for scattering losses through a relaxation time corresponding to , commensurate with experimental data at the considered doping level^{53}. Naturally, resonances are influenced by both the doping level (), the relaxation time (), the dielectric substrate properties (*ε*_{s}), and the characteristic structure dimensions (*L*). For details, we refer to the Supplementary Information.

In the following, we use an efficient method to compute the non-interacting density response matrix , based on Hilbert and fast Fourier transforms (see Ref. ^{40} and Methods section below). Including the effects of a self-consistent Hartree interaction, i.e. within the RPA, the interacting polarizability is given by^{32}

with the Coulomb interaction for *l* ≠ *l*′, and a self-interaction of 0.58 atomic units at *l* = *l*′^{40}. The poles of or equivalently the zeros of the denominator

give the plasmon frequencies. More accurately, since *ε*^{RPA}(*ω*) is a matrix, we in principle seek the eigenvalues *ε _{n}*(

Numerically, the eigenvalues *ε _{n}*(

The calculated eigenvalue loss spectrum for 20 nm graphene equilateral triangles is shown in Figure 1. In the quantum description we distinguish between zigzag and armchair edge terminations, see Supplementary Information. Multiple plasmon peaks are visible in the considered frequency regime. Additionally, at several frequencies, the two considered loss functions (largest and second largest values of ) are nearly identical, while at other frequencies one can be resonant while the other one is not. This is in full accordance with group-theoretical considerations for our structure with *m*-fold rotational symmetry where the *C _{m}* point group leads to either non-degenerate eigenstates or pairs of eigenstates with a double degeneracy

The plasmon modes 1 through 8, being doubly degenerate, are either symmetric or antisymmetric with respect to the mirror symmetry plane. The dipole modes, 1 and 2, with the electric field being polarized orthogonal to each other, are of particular interest due to their strong coupling to optical fields. They can be excited directly by far-field techniques, and the plasmonic local field enhancement is concentrated at the vertices. The modes 3 through 8 penetrate significantly into the bulk, and can be considered as hybridized modes originating from interaction between dipole and bulk modes, because the patterns at the vertices are similar to dipole modes 1 and 2; in addition, the modes 3–6 have finite net dipole momenta, and can couple to far-field radiation. The modes 9–12 are not doubly degenerate, and exhibit threefold rotational symmetry around the center. Although optically dark, these modes are still detectable by suitable near-field techniques. As an example, in an EELS experiment the breathing mode 12 would exhibit the strongest coupling to a nanometer-sized electron beam if this beam passes through the center of the graphene triangle^{58}.

Having described our classical results for graphene triangles, let us now turn to our corresponding tight-binding quantum results. In the quantum description, we calculate the eigenvalue loss spectrum, identify the plasmon mode eigenfrequencies, and then extract the corresponding eigenmodes. Due to the geometrical symmetry, the plasmon eigenmodes should exhibit the same energy degeneracy features as the equilateral triangles in classical calculations, for instance in Figure panels 1(a) and 1(b) several doubly degenerate plasmon modes occur. Figure 3 shows the wave patterns from the quantum calculations, corresponding to the peak labeling in Figure 1(a) and 1(b). We observe that for the armchair case the modes of the same type are blueshifted when compared to their classical counterparts. On the contrary, zigzag termination incur lower plasmon energies with a net redshift compared to the classical case. As an concrete example, the eigenfrequencies of the dipole modes are 0.326 eV, 0.275 eV, and 0.296 eV for the armchair, zigzag, and classical cases, respectively. The associated mode patterns are only slightly different, yet it is clearly seen from the dipole modes, that in zigzag-terminated triangles the mode spreads much more into the bulk while for armchair termination the mode concentrates at the vertices in the same manner as for the classical results. This trend becomes even more evident in the modes 3 and 4 of which the patterns show no hot spots at the vertices. The shifts of armchair and zigzag structures relative to the classical results were recently discussed from an analytical perspective^{39}, and attributed, essentially, to two effects. For the armchair, a nonlocal blueshift accounts for the observed behavior. In the zigzag case, in addition to a nonlocal interaction, the existence of edge states enables an additional dispersive channel, which leads to a net redshift. Similar edge states do not exist for armchair terminations (see Supplementary Figure S3 for additional details). The role of edge states has previously been examined numerically in graphene ribbons^{40}, disks^{41}, and triangles^{43}.

Plasmon hybridization is of both fundamental and practical importance^{59}^{,60}. Hybridization through tuning of the gap distance can be used to achieve better performance through careful design, such as the field enhancement in dimers^{61} and the sensing capabilities in Fano structures^{62}. Here, we study the plasmon hybridization in graphene bow-tie triangles, using the same classical and quantum methods as for individual triangles above. Figure 4 shows the calculated eigenvalue loss spectra for a gap width of 0.5 nm. There are four modes (*n* = 1,2,3,4) in the classical calculations, originating from the four (accounting degeneracy) low-energy dipole modes of the two uncoupled triangles. The hybridization process is illustrated in Figure 5 with a focus on dipole modes, where energies are given with higher precision in order to display the tiny energy shifts associated with the hybridization. We find that each dipole mode in the individual triangles will split into two modes in the bow-tie triangles forming either bonding or antibonding states. The *x*-polarized dipole (0.2964 eV, dipole aligned parallel to bow-tie axis) exhibits large energy splitting, and the corresponding bonding (antisymmetrically coupled) mode has lower energy. However, for the *y*-polarized dipole (0.2963 eV, dipole aligned perpendicular to bow-tie axis) the reduced mode-overlap causes a very small energy splitting. In both cases, the bonding modes are optically active with a net dipole polarization along *x* and *y* direction, respectively.

We find a very similar behavior in the armchair-terminated bow-tie triangles shown in Figure 4a, but with smaller energy splitting, which originates from a weaker mode overlap and weaker coupling strength when compared to the classical calculations. In the zigzag-terminated bow-tie triangles (see Figure 4b), the coupling strength is even weaker and the *x*-polarized dipole exhibits no appreciable energy splitting when compared to the line width of the uncoupled resonances. As a result of this approximate degeneracy, the coupled system exhibits a single broad peak with all four modes merged together. In contrast to the dipole modes, the higher-order plasmon modes show a weak lifting of degeneracy for antisymmetrical and symmetrical states. We mention that the hybridization picture given in Figure 5 is very general, also being satisfied in quantum calculations but with different eigenfrequencies (hybridization diagrams not shown).

The energy splitting or coupling strength depends on the gap width of the bow-tie structures, which can be investigated in the hybridization of *x*-polarized dipoles. We calculate the eigenfrequencies of the hybridized plasmon modes as a function of the width gap, and show the results in Figure 6. The modes in zigzag triangles exhibit very small energy splitting, so we do not show them here. Both in the classical calculations and armchair-terminated triangles, the energy splitting decreases as the gap width increases. The decrease is most pronounced for gap widths below 4 nm, while the variation is weaker for larger separations.

We note that the hybridization of other dimer plasmon modes (other than dipole modes) can be analyzed with a similar result. Generally speaking, the eigenfrequencies of the resulting hybridized modes are decided by two factors: symmetry and coupling strength. Specifically, the antisymmetrically coupled modes (no matter which polarization) have lower energy and modes with less field concentration at the gap region cause weaker coupling and consequently exhibit smaller energy splitting. As a further evidence for this qualitative characterization, we show in Figure 7 the selected twelve plasmon modes from classical calculations, corresponding to the peaks shown in Figure 4c. As compared with Figure 2, they can be understood as linear combinations of the wave patterns in individual structures. Likewise, it is straightforward to envision the wave patterns in armchair and zigzag bow-tie triangles based on the uncoupled modes from Figure 3.

In this article, we have considered and compared classical and quantum aspects of plasmonic eigenmodes in graphene triangular nanostructures. The 2D FEM-approach for calculation of the classical electromagnetic response represents a numerically highly efficient method for electrodynamics in general 2D morphologies of graphene structures in the electrostatic limit (see Supplementary Information for the calculation in hexagonal structures). The simple eigenvalue approach offers a direct pathway to extraction of all plasmonic eigenmodes, not limited to just the optically active, but including also dark modes and highly symmetric breathing modes. The quantum method adopted here is useful for investigating the quantum effects in plasmon excitations of smaller graphene structures, and it offers additional insight into the importance of the particular edge-termination of the underlying atomic lattice. By a sweep of the excitation energy, our calculation of the eigenvalue loss spectra enables direct identification of all plasmonic modes also in the quantum treatment.

We have applied both methods to equilateral triangles, of 20 nm sidelength, both in isolated and in bow-tie configurations. For the isolated nanotriangle we find that the plasmonic response of armchair-terminated triangles is qualitatively similar to the classical case, albeit with a significant and consistent blueshift of all resonances due to nonlocal response. Conversely, the response of zigzag-terminated triangles exhibits several significant differences from its classical counterpart. As a consequence of the existence of localized electronic edge states near zigzag edges, the eigenmodes extend further into the bulk, and are less intense at the vertices. Additionally, we observe a redshift and an pronounced readjustment of the loss-function intensity relative to the classical case.

In the bow-tie configuration we observe plasmon hybridization and associated eigenmode energy splitting, of varying degree depending on treatment; the largest splitting is observed in the classical approach, and the smallest in zigzag structures. Nevertheless, the effects of hybridization are qualitatively similar across the considered cases, with the antisymmetric hybridized modes exhibiting a lowered energy, and with the coupling strength - and associated energy splitting - decreasing when the constituent eigenmodes exhibit lower field intensities in the gap region.

The classical calculations are performed on the two-dimensional domain defined by the geometry of the graphene structure. The domain is discretized using a triangular mesh (see Supplementary Information for details), consisting of a set of elements delimited by a set of vertices . For future reference, we denote the region of the *j*th element by Ω* _{j}*. With a sufficiently dense mesh, a faithful approximation of the Coulomb and Laplacian operators in Eqs. (1) and (2) can then be achieved with the approach described in the following. Crucially, this allows the reduction of the coupled integro-differential equations into simple algebraic equations, as summarized in Eq. (3).

To proceed, we introduce the following notation: the vertices of the *j*th element are denoted as *α _{j}*,

The integration of Eq. (1) can then be approximated directly by a Riemann sum at the centroids. With an aim to ultimately interlink the potentials, * _{k}*, and the densities,

where we have introduced the area of Ω* _{j}* as

Rather than assembling the matrix **B** (the discretized representation of the Laplacian) from Eq. (2) directly, we identify it by its weak form, applying the ideas behind FEM, allowing us to enforce the boundary condition explicitly. Specifically, multiplying onto Eq. (2) an unspecified linear test function , and integrating over the domain, we find

where the boundary condition has been applied at the second equality sign. Next, we specify the linear test function . Concretely, we take the test function as nonzero only on Ω* _{j}*. The value of the test function within Ω

such that equals the approximated () for Ω* _{j}*. Applying Eq. (10) to the left-hand side of Eq. (9) then yields

We can recast this result in terms of the full vertex vectors ** and **** ρ** as equaling

To evaluate the right-hand side of Eq. (9) we require an expression for for , which can be obtained from Eq. (10) using straight-forward algebra, yielding

where denotes a *π*/*2* counterclockwise rotation. Using this, the right-hand side integral in Eq. (9) becomes with **B**^{Rj} denoting a *K* × *K* matrix defined similarly to **B**^{Lj}, i.e. as , but with a *j*-dependent 3 × 3 block matrix

where the subscript-notation *n* + 1 and *n* + 2 indicates forward-cycling by 1 and 2, respectively, in the set {*α _{j}*,

Finally, by summing over all *j*, and defining **B**^{L} Σ_{j}**B**^{Lj} and **B**^{R} Σ_{j}**B**^{Rj}, while noting that the test functions constitute a complete basis in the FEM sense, we can identify the weak form of Eq. (9) as **B**^{L}** ρ** = −

The classical material response is modeled by the bulk conductivity *σ*(*ω*) of graphene through its well-known local-response form^{63}^{,64}

with the first and second terms due to intra- and interband dynamics, respectively. Here, *e* denotes the electron charge, *θ*(*x*) the Heaviside function, and ln(*x*) the natural logarithm.

The tight-binding Hamiltonian for the *π*-electrons is constructed by considering only nearest-neighbor interactions with a hopping strength *t* = 2.8 eV. The associated Hamiltonian matrix-representation is real-valued and symmetric, giving rise to real eigenvalues and eigenvectors.

The direct evaluation of the noninteracting density response matrix of Eq. (4) requires significant computational resources and time, amounting to ~*N*^{4} operations, which must additionally be repeated for each distinct frequency. Significant reduction of computational complexity, to ~*N*^{3}, can be achieved with the aid of Hilbert and fast Fourier transforms (FFT), following a procedure developed in density-functional theory (DFT)^{65}^{,66}, and recently implemented in Ref. ^{40} for the tight-binding model of graphene considered here. We adopt the same technique in our computations.

Furthermore, consideration of the symmetry of , i.e. , leads to an additional reduction of the computational requirements.

We thank Wei Yan, Xiaolong Zhu, and Nicolas Stenger for stimulating discussions. The Center for Nanostructured Graphene is sponsored by the Danish National Research Foundation, Project DNRF58. This work was also supported by the Danish Council for Independent Research–Natural Sciences, Project 1323-00087.

The authors declare no competing financial interests.

**Author Contributions** W.W., T.C. and N.A.M. conceived the basic idea. W.W. performed all numerical simulations. Figures were prepared by W.W. and T.C. All authors interpreted and discussed the results and the writing of the manuscript was done in a joint effort.

- Gramotnev D. K. & Bozhevolnyi S. I. Plasmonics beyond the diffraction limit. Nature Photon. 4, 83–91 (2010).
- Atwater H. A. & Polman A. Plasmonics for improved photovoltaic devices. Nature Mater. 9, 205–213 (2010). [PubMed]
- Barnes W. L., Dereux A. & Ebbesen T. W. Surface plasmon subwavelength optics. Nature 424, 824–830 (2003). [PubMed]
- Ozbay E. Plasmonics: Merging photonics and electronics at nanoscale dimensions. Science 311, 189–193 (2006). [PubMed]
- Oulton R. F., Sorger V. J., Genov D. A., Pile D. F. P. & Zhang X. A hybrid plasmonic waveguide for subwavelength confinement and long-range propagation. Nature Photon. 2, 496–500 (2008).
- Yokogawa S., Burgos S. P. & Atwater H. A. Plasmonic color filters for CMOS image sensor applications. Nano Lett. 12, 4349–4354 (2012). [PubMed]
- Chang W.-S.
*et al.*A plasmonic Fano switch. Nano Lett. 12, 4977–4982 (2012). [PubMed] - Cai W. S., White J. S. & Brongersma M. L. Compact, high-speed and power-efficient electrooptic plasmonic modulators. Nano Lett. 9, 4403–4411 (2009). [PubMed]
- Hess O.
*et al.*Active nanoplasmonic metamaterials. Nature Mater. 11, 573–584 (2012). [PubMed] - West P. R.
*et al.*Searching for better plasmonic materials. Laser Photon. Rev. 4, 795–808 (2010). - Naik G. V. & Boltasseva A. Semiconductors for plasmonics and metamaterials. Phys. Status Solidi-Rapid Res. Lett. 4, 295–297 (2010).
- Tsiatmas A.
*et al.*Superconducting plasmonics and extraordinary transmission. Appl. Phys. Lett. 97, 111106 (2010). - Tassin P., Koschny T., Kafesaki M. & Soukoulis C. M. A comparison of graphene, superconductors and metals as conductors for metamaterials and plasmonics. Nature Photon. 6, 259–264 (2012).
- Ferrari A. C.
*et al.*Science and technology roadmap for graphene, related two-dimensional crystals, and hybrid systems. Nanoscale**7**, 4598-4810 (2015). [PubMed] - Hwang E. H. & Das Sarma S. Dielectric function, screening, and plasmons in two-dimensional graphene. Phys. Rev. B 75, 205418 (2007).
- Jablan M., Buljan H. & Soljačić M. Plasmonics in graphene at infrared frequencies. Phys. Rev. B 80, 245435 (2009).
- Koppens F. H. L., Chang D. E. & García de Abajo F. J. Graphene plasmonics: A platform for strong light-matter interaction. Nano Lett. 11, 3370–3377 (2011). [PubMed]
- Bludov Y. V., Ferreira A., Peres N. M. R. & Vasilevskiy M. I. A primer on surface plasmon-polaritons in graphene. Int. J. Mod. Phys. B 27, 1341001 (2013).
- García de Abajo F. J. Graphene plasmonics: Challenges and opportunities. ACS Photonics 1, 135–152 (2014).
- Zhu X.
*et al.*Experimental observation of plasmons in a graphene monolayer resting on a two-dimensional subwavelength silicon grating. Appl. Phys. Lett. 102, 131101 (2013). - Geim A. K. & Novoselov K. S. The rise of graphene. Nature Mater. 6, 183–191 (2007). [PubMed]
- Castro Neto A. H., Guinea F., Peres N. M. R., Novoselov K. S. & Geim A. K. The electronic properties of graphene. Rev. Mod. Phys. 81, 109–162 (2009).
- Das Sarma S., Adam S., Hwang E. H. & Rossi E. Electronic transport in two-dimensional graphene. Rev. Mod. Phys. 83, 407–470 (2011).
- Yan H.
*et al.*Damping pathways of mid-infrared plasmons in graphene nanostructures. Nature Photon. 7, 394–399 (2013). - Zhu X.
*et al.*Plasmon–phonon coupling in large–area graphene dot and antidot arrays fabricated by nanosphere lithography. Nano Lett. 14, 2907–2913 (2014). [PubMed] - Novoselov K. S.
*et al.*Two-dimensional gas of massless Dirac fermions in graphene. Nature 438, 197–200 (2005). [PubMed] - Zhang Y., Tan Y.-W., Stormer H. L. & Kim P. Experimental observation of the quantum Hall effect and Berry's phase in graphene. Nature 438, 201–204 (2005). [PubMed]
- Ju L.
*et al.*Graphene plasmonics for tunable terahertz metamaterials. Nature Nanotechnol. 6, 630–634 (2011). [PubMed] - Fang Z.
*et al.*Gated tunability and hybridization of localized plasmons in nanostructured graphene. ACS Nano 7, 2388–2395 (2013). [PubMed] - Low T. & Avouris P. Graphene plasmonics for terahertz to mid-infrared applications. ACS Nano 8, 1086–1101 (2014). [PubMed]
- Zhu X., Yan W., Mortensen N. A. & Xiao S. Bends and splitters in graphene nanoribbon waveguides. Opt. Express 21, 3486–3491 (2013). [PubMed]
- Wunsch B., Stauber T., Sols F. & Guinea F. Dynamical polarization of graphene at finite doping. New J. Phys. 8, 318 (2006).
- Christensen J., Manjavacas A., Thongrattanasiri S., Koppens F. H. L. & García de Abajo F. J. Graphene plasmon waveguiding and hybridization in individual and paired nanoribbons. ACS Nano 6, 431–440 (2011). [PubMed]
- Thongrattanasiri S., Koppens F. H. L. & García de Abajo F. J. Complete optical absorption in periodically patterned graphene. Phys. Rev. Lett. 108, 047401 (2012). [PubMed]
- Mortensen N. A., Raza S., Wubs M., Søndergaard T. & Bozhevolnyi S. I. A generalized non-local optical response theory for plasmonic nanostructures. Nature Commun. 5, 3809 (2014). [PubMed]
- Tielrooij K. J.
*et al.*Electrical control of optical emitter relaxation pathways enabled by graphene. Nature Physics**11**, 281-287 (2015). - Christensen T.
*et al.*Nonlocal response of metallic nanospheres probed by light, electrons, and atoms. ACS Nano 8, 1745–1758 (2014). [PubMed] - Wang W. & Kinaret J. M. Plasmons in graphene nanoribbons: Interband transitions and nonlocal effects. Phys. Rev. B 87, 195424 (2013).
- Christensen T., Wang W., Jauho A.-P., Wubs M. & Mortensen N. A. Classical and quantum plasmonics in graphene nanodisks: Role of edge states. Phys. Rev. B 90, 241414(R) (2014).
- Thongrattanasiri S., Manjavacas A. & García de Abajo F. J. Quantum finite-size effects in graphene plasmons. ACS Nano 6, 1766–1775 (2012). [PubMed]
- Thongrattanasiri S. & García de Abajo F. J. Optical field enhancement by strong plasmon interaction in graphene nanostructures. Phys. Rev. Lett. 110, 187401 (2013). [PubMed]
- Manjavacas A. & García de Abajo F. J. Tunable plasmons in atomically thin gold nanodisks. Nature Commun. 5, 3548 (2014). [PubMed]
- Manjavacas A., Thongrattanasiri S. & García de Abajo F. J. Plasmons driven by single electrons in graphene nanoislands. Nanophotonics 2, 139–151 (2013).
- Thongrattanasiri S., Manjavacas A., Nordlander P. & García de Abajo F. J. Quantum junction plasmons in graphene dimers. Laser Photon. Rev. 7, 297–302 (2013).
- Pedersen T. G.
*et al.*Graphene antidot lattices: Designed defects and spin qubits. Phys. Rev. Lett. 100, 136804 (2008). [PubMed] - Nikitin A. Y., Guinea F. & Martín-Moreno L. Resonant plasmonic effects in periodic graphene antidot arrays. Appl. Phys. Lett. 101, 151119 (2012).
- Nikitin A. Y., Guinea F., García-Vidal F. J. & Martín-Moreno L. Edge and waveguide terahertz surface plasmon modes in graphene microribbons. Phys. Rev. B 84, 161407 (2011).
- Nikitin A. Y., Guinea F., García-Vidal F. J. & Martín-Moreno L. Surface plasmon enhanced absorption and suppressed transmission in periodic arrays of graphene ribbons. Phys. Rev. B 85, 081405 (2012).
- Fetter A. L. Magnetoplasmons in a two-dimensional electron fluid: Disk geometry. Phys. Rev. B 33, 5221–5227 (1986). [PubMed]
- Wang W., Apell S. P. & Kinaret J. M. Edge magnetoplasmons and the optical excitations in graphene disks. Phys. Rev. B 86, 125450 (2012).
- Wallace P. R. The band theory of graphite. Phys. Rev. 71, 622–634 (1947).
- Reich S., Maultzsch J., Thomsen C. & Ordejón P. Tight-binding description of graphene. Phys. Rev. B 66, 035412 (2002).
- Tassin P., Koschny T. & Soukoulis C. M. Graphene for terahertz applications. Science 341, 620–621 (2013). [PubMed]
- Andersen K., Jacobsen K. W. & Thygesen K. S. Spatially resolved quantum plasmon modes in metallic nano-films from first-principles. Phys. Rev. B 86, 245129 (2012).
- Andersen K., Jensen K. L., Mortensen N. A. & Thygesen K. S. Visualizing hybridized quantum plasmons in coupled nanowires: From classical to tunneling regime. Phys. Rev. B 87, 235433 (2013).
- Sakoda K.
*Optical Properties of Photonic Crystals*, vol. 80 of*Springer Series in Optical Sciences*(Springer, Berlin, 2005), 2 edn. - Awada C.
*et al.*Selective excitation of plasmon resonances of single Au triangles by polarization-dependent light excitation. J. Phys. Chem. C 116, 14591–14598 (2012). - Schmidt F.-P.
*et al.*Dark plasmonic breathing modes in silver nanodisks. Nano Lett. 12, 5780–5783 (2012). [PMC free article] [PubMed] - Prodan E., Radloff C., Halas N. J. & Nordlander P. A hybridization model for the plasmon response of complex nanostructures. Science 302, 419–422 (2003). [PubMed]
- Nordlander P., Oubre C., Prodan E., Li K. & Stockman M. I. Plasmon hybridization in nanoparticle dimers. Nano Lett. 4, 899–903 (2004).
- Zuloaga J., Prodan E. & Nordlander P. Quantum description of the plasmon resonances of a nanoparticle dimer. Nano Lett. 9, 887–891 (2009). [PubMed]
- Hao F.
*et al.*Symmetry breaking in plasmonic nanocavities: Subradiant LSPR sensing and a tunable Fano resonance. Nano Lett. 8, 3983–3988 (2008). [PubMed] - Hanson G. W. Dyadic Green's functions and guided surface waves for a surface conductivity model of graphene. J. Appl. Phys. 103, 064302 (2008).
- Falkovsky L. A. & Varlamov A. A. Space-time dispersion of graphene conductivity. Eur. Phys. J. B 56, 281–284 (2007).
- Shishkin M. & Kresse G. Implementation and performance of the frequency-dependent GW method within the PAW framework. Phys. Rev. B 74, 035101 (2006).
- Miyake T. & Aryasetiawan F. Efficient algorithm for calculating noninteracting frequency-dependent linear response functions. Phys. Rev. B 61, 7172–7175 (2000).

Articles from Scientific Reports are provided here courtesy of **Nature Publishing Group**