PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
 
Sci Rep. 2015; 5: 9535.
Published online 2015 April 9. doi:  10.1038/srep09535
PMCID: PMC5396075

Plasmonic eigenmodes in individual and bow-tie graphene nanotriangles

Abstract

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 limit1, squeezing the light down to nanoscale2, and large local field enhancement3. 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 scales4. Developments in nanofabrication technology have stimulated a series of plasmon-based devices like waveguides5, filters6, switches7, and modulators8. 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 plasmonics9. Alternatively, attempts have already been made to study plasmonics in materials other than metals10, for example doped semiconductors11 and superconductors12,13.

Graphene and other low-dimensional crystals are now emerging as interesting materials for exciting science and technology14. 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 well15,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 approaches20. The charge carriers in graphene obey linear energy dispersion at lower energies close to the Dirac points, thus resembling the linear dispersion of photons21,22,23. Experimental investigations of carrier transport show that the mobility limited by impurity scattering can exceed 15.000 cm2/Vs at room temperature21, which gives the intrinsic loss in graphene one order of magnitude less than the noble metals. Despite relaxation due to phonon scattering24,25, graphene achieves superior plasmonic performance in propagation length and field enhancement26,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 experimentally28,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 technology30. 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 loss31.

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 representation33,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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m14.jpg = 0.4 eV), the common assumption An external file that holds a picture, illustration, etc.
Object name is srep09535-m15.jpg 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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m16.jpg holds, while instead the tiny gap (An external file that holds a picture, illustration, etc.
Object name is srep09535-m17.jpg) promotes nonlocal effects35. For optical excitation in the near field17,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 scale37. In this regime, both semiclassical hydrodynamic38,39 and full quantum approaches have been proposed40,41, similar to those recently developed for metals35,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 regime43. We will emphasize dimers formed by electrically mutually disconnected graphene islands, while graphene dimers connected by quantum junctions44 and extended complimentary structures (electronic45 and plasmonic25,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 channels39, though we will not pursue this scheme for our classical calculations presently.

In our quantum treatment, we employ a tight-binding description40,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 nanoparticles35.

Results

Classical Description

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 tg. This assumed 3D film has an effective bulk permittivity, ε(ω) = ε0 + (ω)/(ωtg), where ε0 is the vacuum permability and σ(ω) denotes the surface conductivity as obtained from e.g. the local-response limit of the RPA33,34. Evidently, the artificial thickness tg 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 structure17,36) would necessarily converge to physically meaningful quantities in the limit tg → 0. Alternatively, in nanostructures with high symmetry, e.g. in ribbons47,48 or disks39,49, one may take advantage of modal expansion methods39,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 [var phi] and the induced density ρ. In particular, we note that the total potential [var phi](An external file that holds a picture, illustration, etc.
Object name is srep09535-m110.jpg) is governed by Coulomb's law

An external file that holds a picture, illustration, etc.
Object name is srep09535-m1.jpg

where [var phi]ext(An external file that holds a picture, illustration, etc.
Object name is srep09535-m111.jpg) denotes the external potential, L is an auxiliary quantity such as the feature length of the structure which makes the surface integral dimensionless, ρ(An external file that holds a picture, illustration, etc.
Object name is srep09535-m112.jpg′) 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 J2D = −σ(ω)[nabla]2D[var phi](An external file that holds a picture, illustration, etc.
Object name is srep09535-m114.jpg) into the continuity equation iωρ(An external file that holds a picture, illustration, etc.
Object name is srep09535-m113.jpg) = [nabla]2D·J2D(An external file that holds a picture, illustration, etc.
Object name is srep09535-m119.jpg), which for r restricted to the plane of the graphene structure gives

An external file that holds a picture, illustration, etc.
Object name is srep09535-m2.jpg

with An external file that holds a picture, illustration, etc.
Object name is srep09535-m18.jpg the 2D Laplace operator. Equation (2) is solved subject to the assumption of charge neutrality, i.e. An external file that holds a picture, illustration, etc.
Object name is srep09535-m19.jpg, implying that An external file that holds a picture, illustration, etc.
Object name is srep09535-m20.jpg on the boundary of the domain, with An external file that holds a picture, illustration, etc.
Object name is srep09535-m21.jpg 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 [var phi] 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 [var phi] = [var phi]ext + (4πεsL)−1Aρ and ρ = (ω)ω−1B[var phi], which we combine to get

An external file that holds a picture, illustration, etc.
Object name is srep09535-m3.jpg

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(ω) = (ω)/(4πεs) is a geometry-independent scalar29. The term in the square brackets on the left-hand side of Eq. (3) represents the matrix-form of the effective (classical) frequency-dependent dielectric function εCLA(ω) [equivalent] 1f(ω)BA, connecting total and external potentials via εCLA(ω)[var phi] = [var phi]ext. In the absence of an external potential ([var phi]ext = 0), Eq. (3) becomes an eigenvalue problem for the matrix BA. The resulting eigenvalues λn are associated with plasmon frequencies ωn through An external file that holds a picture, illustration, etc.
Object name is srep09535-m22.jpg, and the associated eigenvectors are induced charge densities ρn in a finite-element representation. The corresponding eigenpotentials are denoted as [var phi]n, and within the graphene plane they can be computed directly as [var phi]n = Aρn. Following this classical approach, all plasmonic eigenmodes for a specific structure can be obtained as the solution of a single eigenvalue problem. This constitutes an attractive computational approach that can give direct insight in the plasmonic eigenstates that one would be able to probe with various experimental techniques.

Quantum Mechanical Tight-Binding Description

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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m23.jpg orbitals are well separated in energy from the four σ bands arising from sp2 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 approximation51,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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m24.jpg 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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m25.jpg, is then built from the electronic states whose elements are given by15,16,32

An external file that holds a picture, illustration, etc.
Object name is srep09535-m4.jpg

where An external file that holds a picture, illustration, etc.
Object name is srep09535-m26.jpg denotes the Fermi–Dirac distribution function associated with the state with energy An external file that holds a picture, illustration, etc.
Object name is srep09535-m27.jpg and wave function An external file that holds a picture, illustration, etc.
Object name is srep09535-m28.jpg (l labels each of the carbon atoms), while kB and An external file that holds a picture, illustration, etc.
Object name is srep09535-m29.jpg 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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m30.jpg and a temperature T = 300 K corresponding to a thermal energy of An external file that holds a picture, illustration, etc.
Object name is srep09535-m31.jpg. We phenomenologically account for scattering losses through a relaxation time An external file that holds a picture, illustration, etc.
Object name is srep09535-m32.jpg corresponding to An external file that holds a picture, illustration, etc.
Object name is srep09535-m33.jpg, commensurate with experimental data at the considered doping level53. Naturally, resonances are influenced by both the doping level (An external file that holds a picture, illustration, etc.
Object name is srep09535-m34.jpg), the relaxation time (An external file that holds a picture, illustration, etc.
Object name is srep09535-m35.jpg), 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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m36.jpg, 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 by32

An external file that holds a picture, illustration, etc.
Object name is srep09535-m5.jpg

with the Coulomb interaction An external file that holds a picture, illustration, etc.
Object name is srep09535-m37.jpg for ll′, and a self-interaction of 0.58 atomic units at l = l40. The poles of An external file that holds a picture, illustration, etc.
Object name is srep09535-m38.jpg or equivalently the zeros of the denominator

An external file that holds a picture, illustration, etc.
Object name is srep09535-m6.jpg

give the plasmon frequencies. More accurately, since εRPA(ω) is a matrix, we in principle seek the eigenvalues εn(ω) of the εRPA(ω) whose real part approaches zero54. In practice there is also loss, for example due to An external file that holds a picture, illustration, etc.
Object name is srep09535-m39.jpg in An external file that holds a picture, illustration, etc.
Object name is srep09535-m40.jpg. The eigenfrequencies εn(ω) are therefore complex-valued, with the imaginary parts denoting the plasmon peak broadening. With this in mind, we finally define the plasmon resonance frequencies from the local maxima of An external file that holds a picture, illustration, etc.
Object name is srep09535-m41.jpg54,55.

Numerically, the eigenvalues εn(ω) are obtained by diagonalizing the RPA dielectric function εRPA(ω) for each frequency. An N-atom nanostructure entails n = 1,2,3, …, N distinct eigenvalues. Out of these we focus in the following on eigenvalues with largest and second-largest values of An external file that holds a picture, illustration, etc.
Object name is srep09535-m42.jpg. Their corresponding eigenvectors are the induced charge densities ρn, and similarly the eigenpotentials [var phi]n can be obtained by performing the Coulomb integral. For comparison with the quantum treatment, we also calculate the eigenvalue loss spectrum in the classical framework by carrying out diagonalization of the classical effective dielectric function εCLA(ω).

Plasmonic Eigenmodes in Individual Triangles

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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m43.jpg) 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 Cm point group leads to either non-degenerate eigenstates or pairs of eigenstates with a double degeneracy56. The degeneracy can be explored further by considering the eigenmodes, expressed e.g. by the in-plane potential, and in particular their symmetries. In the classical approach, the eigenmodes appear as eigenvectors of the matrix BA of Eq. (3). Considering the two lowest eigenstates causing the resonance around 0.3 eV in Figure 1(c), we numerically find the eigenfrequencies to be 0.2964 eV and 0.2963 eV. The small energy difference of 0.1 meV illustrates the numerical accuracy (symmetry breaking) associated with the fact that our finite-element mesh does not comply with the threefold rotational symmetry of the graphene triangle. In Figure 2 we show corresponding in-plane potential distributions of the twelve lowest-energy eigenmodes, again calculated in the classical framework. The eigenmodes are responsible for the primary features of Figure 1c; specifically, the loss-function exhibits peaks at the resonance energies of the eigenmodes. The peaks are each assigned a label (n = 1,2,3,…), corresponding to the eigenmode enumeration in Figure 2. A one-to-one correspondence is evident and whenever the spectrum in Figure 1 suggests a pair of degenerate states, the corresponding modes in Figure 2 support that they are indeed pairs of orthogonal and degenerate states. The energy degeneracies exhibited here are a direct consequence of the symmetries of the considered nanostructure, as required by group theory57.

Figure 1
Eigenvalue loss spectrum in graphene triangles.
Figure 2
In-plane classical eigenmode potentials.

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 triangle58.

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 perspective39, 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 ribbons40, disks41, and triangles43.

Figure 3
In-plane quantum eigenmode potentials.

Plasmon Hybridization in Bow-Tie Triangles

Plasmon hybridization is of both fundamental and practical importance59,60. Hybridization through tuning of the gap distance can be used to achieve better performance through careful design, such as the field enhancement in dimers61 and the sensing capabilities in Fano structures62. 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.

Figure 4
Eigenvalue loss spectrum in graphene dimers.
Figure 5
Hybridization of plasmons in graphene dimers.

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.

Figure 6
Gap dependence of hybridization in graphene dimers.

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.

Figure 7
In-plane classical eigenmode potentials of graphene dimers.

Discussion

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.

Methods

Classical Calculations

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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m44.jpg delimited by a set of vertices An external file that holds a picture, illustration, etc.
Object name is srep09535-m45.jpg. For future reference, we denote the region of the jth 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 jth element are denoted as αj, βj, and δj, and together define the element centroid coordinates, An external file that holds a picture, illustration, etc.
Object name is srep09535-m46.jpg and values An external file that holds a picture, illustration, etc.
Object name is srep09535-m47.jpg (representing e.g. the density)

An external file that holds a picture, illustration, etc.
Object name is srep09535-m7.jpg

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, [var phi]k, and the densities, ρk, at all vertex positions An external file that holds a picture, illustration, etc.
Object name is srep09535-m48.jpg, we can then evaluate Eq. (1) at the kth vertex to find

An external file that holds a picture, illustration, etc.
Object name is srep09535-m8.jpg

where we have introduced the area of Ωj as sj. Note that no explicit handling of a semi-divergent self-interaction is necessary since An external file that holds a picture, illustration, etc.
Object name is srep09535-m49.jpg for all j. The matrix A can then be assembled by noting its equivalent definition in element form, An external file that holds a picture, illustration, etc.
Object name is srep09535-m50.jpg, which upon comparison with Eq. (8) allows the identification An external file that holds a picture, illustration, etc.
Object name is srep09535-m51.jpg, with the m-sum running over the vertex sites of Ωj, i.e. m [set membership] {αj,βj,δj}, and where the twice-subscripted δ denotes the Kronecker delta.

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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m52.jpg explicitly. Specifically, multiplying onto Eq. (2) an unspecified linear test function An external file that holds a picture, illustration, etc.
Object name is srep09535-m53.jpg, and integrating over the domain, we find

An external file that holds a picture, illustration, etc.
Object name is srep09535-m9.jpg

where the boundary condition An external file that holds a picture, illustration, etc.
Object name is srep09535-m54.jpg has been applied at the second equality sign. Next, we specify the linear test function An external file that holds a picture, illustration, etc.
Object name is srep09535-m55.jpg. Concretely, we take the test function as nonzero only on Ωj. The value of the test function within Ωj is most conveniently specified in a local barycentric coordinate system with parameters η [set membership] [0,1] and An external file that holds a picture, illustration, etc.
Object name is srep09535-m56.jpg, within which An external file that holds a picture, illustration, etc.
Object name is srep09535-m57.jpg for An external file that holds a picture, illustration, etc.
Object name is srep09535-m116.jpg [set membership] Ωj. Within each region Ωj we interpolate function values An external file that holds a picture, illustration, etc.
Object name is srep09535-m58.jpg, e.g. potentials or densities, according to their associated vertex values

An external file that holds a picture, illustration, etc.
Object name is srep09535-m10.jpg

such that An external file that holds a picture, illustration, etc.
Object name is srep09535-m59.jpg equals the approximated [var phi](An external file that holds a picture, illustration, etc.
Object name is srep09535-m115.jpg) for An external file that holds a picture, illustration, etc.
Object name is srep09535-m117.jpg [set membership] Ωj. Applying Eq. (10) to the left-hand side of Eq. (9) then yields

An external file that holds a picture, illustration, etc.
Object name is srep09535-m60.jpg

We can recast this result in terms of the full vertex vectors [var phi] and ρ as equaling [var phi]TBLjρ, where BLj denotes a K × K matrix defined by An external file that holds a picture, illustration, etc.
Object name is srep09535-m61.jpg with the nm-sum restricted to {n,m} [set membership] {αj, βj, δj}, with nonzero elements defined by a 3 × 3 block matrix An external file that holds a picture, illustration, etc.
Object name is srep09535-m62.jpg.

To evaluate the right-hand side of Eq. (9) we require an expression for An external file that holds a picture, illustration, etc.
Object name is srep09535-m63.jpg for An external file that holds a picture, illustration, etc.
Object name is srep09535-m118.jpg, which can be obtained from Eq. (10) using straight-forward algebra, yielding

An external file that holds a picture, illustration, etc.
Object name is srep09535-m11.jpg

where An external file that holds a picture, illustration, etc.
Object name is srep09535-m64.jpg denotes a π/2 counterclockwise rotation. Using this, the right-hand side integral in Eq. (9) becomes An external file that holds a picture, illustration, etc.
Object name is srep09535-m65.jpg with BRj denoting a K × K matrix defined similarly to BLj, i.e. as An external file that holds a picture, illustration, etc.
Object name is srep09535-m66.jpg, but with a j-dependent 3 × 3 block matrix

An external file that holds a picture, illustration, etc.
Object name is srep09535-m12.jpg

where the subscript-notation n + 1 and n + 2 indicates forward-cycling by 1 and 2, respectively, in the set {αj, βj, δj} (e.g. if n = βj then n + 1 = δj and n + 2 = αj, and equivalently for m).

Finally, by summing over all j, and defining BL [equivalent] ΣjBLj and BR [equivalent] ΣjBRj, while noting that the test functions An external file that holds a picture, illustration, etc.
Object name is srep09535-m67.jpg constitute a complete basis in the FEM sense, we can identify the weak form of Eq. (9) as BLρ = −(ω)ω−1BR[var phi], from which we identify the desired matrix B with the property ρ = (ω)ω−1B[var phi] as B = −(BL)−1BR.

The classical material response is modeled by the bulk conductivity σ(ω) of graphene through its well-known local-response form63,64

An external file that holds a picture, illustration, etc.
Object name is srep09535-m13.jpg

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.

Quantum Calculations

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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m68.jpg of Eq. (4) requires significant computational resources and time, amounting to ~N4 operations, which must additionally be repeated for each distinct frequency. Significant reduction of computational complexity, to ~N3, 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 An external file that holds a picture, illustration, etc.
Object name is srep09535-m69.jpg, i.e. An external file that holds a picture, illustration, etc.
Object name is srep09535-m70.jpg, leads to an additional reduction of the computational requirements.

Supplementary Material

Supplementary Information:

Supplementary Information

Acknowledgments

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.

Footnotes

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.

References

  • 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