Highly resonant graphene plasmon hotspots in complex nanoresonator geometries

Van der Waals surface polariton nanostructures are promising candidates for miniaturisation of electromagnetic devices through the nanoscale confinement of infrared light. To fully exploit these nanoresonators, a computationally efficient model is necessary to predict polariton behaviour in complex geometries. Here, we develop a general wave model of surface polaritons in 2D geometries smaller than the polariton wavelength. Using geometric approximation widely tuneable infrared nanoimaging and local work function microscopy, we test this model against complex mono-/bi-layer graphene plasmon nanoresonators. Direct imaging of highly resonant graphene plasmon hotspots confirms that the model provides quantitatively accurate, analytical predictions of nanoresonator behaviour. The insights built with such models are crucial to the development of practical plasmonic nanodevices.

Recently, scattering scanning near-field optical microscopy (s-SNOM) has enabled direct imaging of plasmon modes in nanoresonators [17,19]. Until now, the theoretical understanding of such experiments has required time-consuming and computationally expensive simulations using finite or boundary element methods [16,17,19], which require the input of a free parameter, the Fermi energy. These drawbacks mean that only simple 2D geometries, such as discs and rectangles [16,19] or ribbons [11,17,22], have been studied. Furthermore, the reliance on numerical simulations inhibits an intuitive understanding of polariton nanoresonators.
Here, we introduce a simple wave model of how polaritons, including plasmons, behave in 2D geometries of length scales, L, smaller than the polariton wavelength (i.e. L < λ sp ). This model allows simple analytical calculations, which provide quantitatively accurate predictions of polariton behaviour in nanoresonators.
As the model requires prior knowledge of the plasmon dispersion relation, we first determine exper imentally the Fermi energy of the single-layer graphene (1LG) in the sample under study. Typical methods for measuring E F are the Hall effect, Kelvin probe force microscopy (KPFM), vector decomposition of Raman spectra, ultraviolet photoemission spectr oscopy (UPS) and angle-resolved photoemission spectr oscopy (ARPES). However, each method has its own set of advantages and drawbacks. For example, Hall effect measurements offer straightforward and direct measurement of the carrier density, which can then be used to calculate the Fermi energy. However, the Hall effect method typically requires fabrication of electrical devices with the typical sizes ranging from hundreds of nanometres to several millimeters [23,24]. This leads to averaging of the Hall voltage measurement over the large sample area and can therefore give a misrepresented E F for non-uniform graphene samples.
KPFM is often used to map the local work function of graphene, which in turn can be used to estimate E F . However, KPFM is extremely sensitive to surface contamination as well as substrate and environmental doping [25,26], and estimation of the sample's work function requires meticulous calibration of the KPFM probe [27].
The vector decomposition of the Raman G and 2D modes can also be used to determine the carrier density of graphene, but this method is currently limited to p-type 1LG samples [28]. Moreover, the spatial resolution of Raman mapping is limited by the relatively large spot size of the laser (100s of nanometres), making it unsuitable for samples with non-uniform layer thicknesses.
Finally, nano-ARPES offers direct measurement of E F [29], but the technique requires ultrahigh vacuum conditions that significantly change the doping of graphene, thus making it difficult to form direct compariso ns to ambient conditions. In 2012, Chen et al [11] demonstrated that the wavelength of graphene plasmons can be measured by analysing the surface plasmon interference fringes that appear in images obtained by scattering-type scanning near-field optical microscopy (s-SNOM). Using a first order approximation of the plasmon dispersion relation, they inferred an approximate value for the Fermi energy of exfoliated graphene, E F ≈ 400 meV. Since then, many experiments have used this technique to study graphene and other layered 2D materials with applications in the study of chemical doping [30], grain boundary phenomena [31], graphene plasmonic nano-resonators [19,32], and hybrid systems with hyperbolic [33] or phonon polaritons [34]. However, the accuracy and energy range of these studies have so far been limited to a small number of discrete lines in a narrow spectral range, corresponding to freespace wavelengths of λ 0 ≈ 9.2 − 10.6 µm that can be accessed by the CO 2 lasers typically used for s-SNOM measurements [11,35].
Here, we demonstrate broad measurement of plasmon dispersion with high spectral resolution to provide a high-accuracy, quasi-local measurement of the graphene Fermi energy. This enables accurate prediction of plasmon hotspots in mono-/bi-layer graphene constrictions using our model. We confirm these predictions using s-SNOM to directly image the behaviour of plasmons in arbitrary sub-wavelength (L < λ sp ) constrictions. These new experimentally confirmed, analytical insights into the behaviour of nano-constrained 2D polariton hotspots open the path to development of more complex nanoresonator devices.

Plasmon nano-imaging
S-SNOM is a type of scanning probe microscopy that relies on 'tapping mode' atomic force microscopy (AFM), i.e. with a vertical dither applied to the probe at its mechanical resonance frequency ( f 0 ≈ 280 kHz), which tracks the sample surface height. When the probe is in a close vicinity of the sample, the sharp platinum-coated tip radially scatters incident laser light into the graphene surface plasmon modes, which then propagate in the sample plane (figure 1). When the plasmons encounter a defect, for example the edge of the 2D graphene layer, they are partially reflected back and produce a standing wave interference pattern, whose fringes have a spatial periodicity given by half the surface plasmon wavelength, λ sp /2 [11]. The amount of energy subsequently scattered out of the plasmon modes and back into the s-SNOM detection system is proportional to the field intensity at the probe, allowing the fringes to be mapped out in the s-SNOM image. The backscattered signal is demodulated at a higher harmonic (n > 1) of f 0 . This eliminates signals due to background scattered light, and makes the measurement sensitive only to light scattered by its interaction with the near-field intensity in the nanoscale region between the probe apex and the sample [36]. All s-SNOM measurements presented in this work use the near-field scattered amplitude at the second harmonic, s 2 .
Here, we pair s-SNOM with a radiation source consisting of a bank of widely tuneable quantum cascade lasers (QCL; MIRcat, Daylight Solutions). This allows us the capability to tune through a free-space wavelength range λ 0 = 8.93-10.43 µm with high spectral resolution. In turn, this enables the full graphene plasmon dispersion curve to be obtained and highly resonant plasmonic phenomena investigated.

Sample preparation
The graphene used in the present study was grown on semi-insulating 4H-SiC(0 0 0 1) substrates in a hot-wall reactor (Aixtron VP508). First, a complete interfacial layer (IFL) was grown under an argon laminar flow at 1600 °C. Following IFL growth, the sample was annealed in hydrogen-rich environment. The hydrogen penetrates between the IFL and substrate, forming Si-H bonds, transforming the IFL to quasi-free-standing (QFS) 1LG. For more details, see [29]. As it is not possible to entirely avoid nucleation growth in certain areas, e.g. near the edges of SiC terraces, such areas are covered with bi-layer graphene (2LG). QFS graphene on SiC produced with this growth method is typically p-type doped, which is attributed to the spontaneous polarisation of the hexagonal SiC [37,38].
To investigate the electronic properties of the sample sheet, a van der Pauw structure was fabricated. Firstly, using electron beam lithography, the sample was etched into a 10 µm × 10 µm square, ensuring the external regions were etched down to the SiC substrate. Four gold contacts were then deposited onto the corners of the graphene using electron beam evaporation. The contacts were annealed in argon at 400 °C for 10 min to improve adhesion of the contacts to the graphene. Finally, large pads were fabricated using a lift-off process to allow for electrical connections. The sample preparation process is detailed in [39]. The entire graphene device was cleaned mechanically of resist residues using contact-mode AFM. Figure 2(a) shows the topography of the sample obtained with AFM, featuring height steps (2-6 nm) from the underlying SiC substrate, together with particulate remnants from the sample fabrication processes (small white spots typically seen on the SiC substrate). Figure 2(b) shows the work function, Φ, of the same area. This was obtained by first calibrating the KPFM probe's work function, Φ Probe , against gold, Φ Au = 4.82 eV, using Φ Probe = Φ Au + eV SP (Au), where V SP is the surface potential measured on gold, followed by applying Φ Sample = Φ Probe − eV SP (Sample) to the surface potential image to obtain the work function of the sample [27]. Given that KPFM is sensitive only to changes in the surface potential, it clearly distinguishes the structure of 1LG, 2LG, SiC and Au due to their inherent differences in work function. For 1LG and 2LG, the work functions are Φ 1LG = 4.81 eV and Φ 2LG = 4.56 eV, respectively, and by assuming the intrinsic work function of graphene as Φ 0 = 4.47 ± 0.05 eV from our previous work [26], we estimate the Fermi energy for 1LG as E KPFM F = 340 ± 10 meV. The Fermi energy can also be measured via the Hall effect, as E 1LG

Measurement of Fermi energy with KPFM and the Hall effect
√ πn e for single-layer graphene (1LG), where ν F is the Fermi velocity, is reduced Planck's constant and n e is the carrier density of graphene. For bi-layer graphene (2LG), the Fermi energy is given by E 2LG F = π 2 n e /2m * , where m * is the effective mass of the charge carriers. Using this method, we measure the Fermi energy (in the van der Pauw geometry [39]) to be E Hall F = 333 ± 1 meV, which compares well with the KPFM measurement above. The assignment of 2LG in the middle of the graphene sheet was also confirmed by mapping the Raman 2D-peak width, where it was determined as ~60 cm −1 for 2LG versus ~20 cm −1 for 1LG (figure 2(e)) [40]. Furthermore, the 2D peak can be fitted by a single Lorentzian for 1LG and a broad peak with shoulders requiring four Lorentzians, which is characteristic for 2LG (figure 2(f)) [40][41][42].

High-accuracy measurement of Fermi energy using plasmon nano-imaging
The s-SNOM amplitude, s 2 , image in figure 2(c) shows a strong mid-IR signal from the 1LG and reveals plasmon interference fringes at the 1LG/SiC and 1LG/2LG boundaries. The s-SNOM amplitude image shows interference fringes emanating from the right and top of the image (as indicated by the white arrows) that arise from the strong plasmon reflection from the straight, lithographically defined edges of the graphene sheet (i.e. the 1LG/SiC boundary). For homogeneous 1LG, the full dispersion relation of plasmons as a function of optical frequency, ω, is given by [1]: where ε sub is the substrate dielectric function, k 0 ≡ 2π λ0 = ω/c is the free space wavevector, k sp ≡ 2π/λ sp is the in-plane plasmon wavevector and σ(ω) is the graphene conductivity. Under the random phase approximation, and for excitation energies sufficiently below E F so as to avoid interband transitions ( ω 2E F ), the graphene conductivity depends linearly on the Fermi energy E F according to [1]: where e is the electronic charge and τ is the relaxation time of charge carriers in the graphene. By smoothly tuning the QCL through a wide range from λ 0 = 8.93 to 10.43 µm, in steps of 100 nm, we map out the variation of the fringe spacing with high spectral resolution. Figures 3(a) and (b) show maps of the plasmon interference fringes at λ 0 = 10.03 and 9.83 µm, respectively. Obtaining maps throughout the specified range of the QCL, we determine the graphene plasmon dispersion over more than a decade in the surface plasmon wavelength, λ sp = 140-1700 nm ( figure 3(c)). As expected, the fringe spacing, and so the plasmon wavelength, increases with increasing λ 0 . Despite the larger error bar for the final data point, which stems from the poorer contrast of the peaks and troughs of the plasmon standing wave, we can extract the precise value of E F by combining equations (1) and (2) and fitting the experimental data with an errorweighted numerical minimisation. Across the spectral range used here (λ 0 = 8.93 to 10.43 µm), the dielectric function of the SiC substrate decreases monotonically from ε sub = 3.3 to −0.5, and the exact form is taken from the literature [11]. The resulting fitted dispersion relation yields E SNOM F = 298 ± 4 meV ( figure 3(c)).
The larger estimation of E F from the Hall effect is attributed to the averaging from the presence of a small quantity of 2LG, which typically exhibits a higher carrier density than 1LG [40]. Furthermore, the Hall effect and KPFM measurements were performed immediately after mechanical cleaning, whereas s-SNOM measurements were made after the sample had chemically stabilised to ambient conditions. The observed variations between E Hall F and E KPFM F , and E SNOM F are consistent with the variation of the carrier concentration and work function of graphene due to prolonged exposure to the environment [43].

Comparison of plasmon reflection at different boundaries
Additional fringes are also present from plasmons scattering at the boundaries between 1LG and 2LG (see figure 4(a), obtained at λ 0 = 9.68 µm). These reflections are generally weaker than those at the edges of the 1LG sheets (i.e. from 1LG/SiC boundary). Chen et al reported on plasmons in graphene reflecting purely from the terrace step edges in SiC [35]. In our case, although there is an additional complication of 2LG also being present at step edges, we believe the latter is not the primary cause of plasmon reflection at 1LG/2LG boundaries, given that we also observe plasmon reflections from smaller 2LG islands (blue arrows in figure 2(c)), where the step height is ~0.3 nm ( figure 2(d)). The decrease in surface height from 1LG to the islands is consistent with epitaxial growth of 2LG rather than an exposed region of SiC, which would show as higher [39]. Furthermore, as discussed in the supplementary material (stacks.iop.org/ TDM/6/021003/mmedia), the strengths and shapes of the interference fringes around these islands are in excellent agreement with the fringes associated with the Raman-confirmed 1LG/2LG boundaries (figure 2(e)), and inconsistent with those at the 1LG/SiC boundary. In this case, the plasmon reflection occurs due to the change in carrier density across the 1LG/2LG boundary ( figure 4(b)), as revealed by Fei et al, where they apply a back gate voltage to vary the Fermi energy of single-, bi-and tri-layer graphene [39,[44][45][46].
Finally, we observe a multitude of plasmon focusing hotspots (indicated by green arrows) in 1LG constrained by the 2LG 'strips' in the middle of figure 4(a). These strips of 2LG are known to form, as in this instance, near to terraces in the SiC substrate, as illustrated in figure 4(c). These hotspots arise from the extreme spatial confinement of plasmons in structures smaller than the plasmon wavelength, λ sp .

Plasmon hotspots in sub-wavelength geometries 3.2.1. Theoretical model of surface wave hotspots in sub-wavelength 2D geometries
To understand the nature of the plasmon hotspots that appear in figure 4(a), we propose a theoretical model of  how plasmons, and indeed any surface wave including phonon-polaritons, behave in constrictions of length scales, L, smaller than the plasmon wavelength (i.e. L < λ sp ). To that end, figure 5(a) provides definitions of salient features of fringes caused by plasmon reflection from an arbitrary boundary, wherein the nth fringe has a strength, F n , above a background near-field amplitude, s 2,bg . The first fringe is located at a distance L 0 from the boundary, and the spacing between fringes n and n + 1 is some distance L n . According to the literature [11,34], and consistent with our own experimental measurements above, all fringe-to-fringe spacings that arise from the interference of tip-launched and reflected plasmons are equal to half the plasmon wavelength, i.e. L n = λ sp /2 ∀ n 1. Furthermore, fringe heights typically decay exponentially in a lossy medium [30]. Therefore, the strengths of successive fringes decrease in a geometric progression, and the ratio of fringe strengths, ρ n ≡ F n+1 /F n , is well approximated by the relative strengths of the first two fringes, i.e. that ρ n ≈ ρ ≡ F 2 /F 1 .
First, we consider the simplest sub-wavelength geometry, that of a cavity comprising two parallel opposing boundaries, as illustrated in figure 5(b). As the first fringe is typically much stronger than subsequent fringes (i.e. ρ 1), a strong central hotspot will form at a critical cavity width, L c, = 2L 0 , as this ensures that the fringes from plasmons reflected at either boundary will overlap exactly, in a 2D analogue of a Fabry-Pérot resonator. Any additional contribution from higher order fringes (n > 1) will depend on the size of L 0 relative to L n 1 . For these higher order fringes to constructively interfere with the central hotspot, it is required that L 0 = m ). This means that only the (mn + 1) th fringes constructively interfere with the first fringe to form a central hotspot with strength where ρ < 1. Figure 5(c) illustrates an extension of this analysis to 2D cavities with regular polygonal geometries. As in the previous example, a strong central hotspot will form when first fringe of constructive plasmon interference from each boundary coincides with the geometrical centre of the polygon. As such, the critical side length, L c,p , of a p-sided polygon can be determined from simple trigonometric considerations. For example, an equilateral triangle has a critical side length of L c,3 = 2 √ 3L 0 , and for a square L c,4 = 2L 0 . Generally, a regular polygon with p sides has L c,p = 2L 0 tan (π/p ) .
In the limit of a circle (p → ∞), equation (4) yields L c,∞ = 0, as expected since a circle has zero side length. Therefore, it is necessary to redefine L c,∞ as the critical diameter of a circle of radius L 0 , so that L c,∞ = 2L 0 . In general, the central hotspots in such 2D constrictions will have greater strengths than that of a 1D cavity discussed above. The final geometry we consider here is a wedge with angle θ between the two boundaries, as shown in figure 5(d). This geometry will create a hotspot at the intersection of the first fringe due to each boundary. This hotspot will occur at a distance from the point of the wedge. This corresponds to a critical wedge width as above) exhibits a strong hotspot at wedge width W ∼ 0.6λ sp [11]. This is in good agreement with the theoretical relationship W ≈ 0.52λ sp predicted by equation (6). Using these simple geometric building blocks, it should be straightforward to extend this model to more complicated composite geometries.

Experimental validation of hotspot model in graphene
The magnified KPFM image in figure 6(a) shows the 2LG strips in greater detail and identifies four geometric shapes (numbered 1-4) that will be used to test the model introduced in the previous section.
The KPFM is critical for determining the geometries of these shapes, as the 1LG/2LG boundaries are barely visible in the topography. Figures 6(b)-(f) show plasmon nanoimaging of this region across a range of incident wavelengths from λ 0 = 9.78 µm to 9.38 µm using the fine tuneability of our setup. Figure 6(g) shows that the hotspot in region 1 moves towards the wide end of the θ ≈ 30 • wedge (see figure 5(d)) as λ 0 increases. Using equation (5) and L 0 = λsp 2 from figure 2 gives theoretical hotspot position y ≈ 1.9λ sp , which is in excellent agreement with the experimental relationship L 0 = (1.6 ± 0.4) λ sp taken from the linear fit shown in figure 5(k).
Region 2 can be modelled approximately as a 1D cavity ( figure 5(b)) of width L c = 200 ± 50 nm (from KPFM measurements in the supplementary material) in addition to connecting 1LG/2LG boundary. KPFM is required to map the geometry of our constrictions, as there is too little contrast between 1LG and 2LG in the topography. Such a structure is predicted to have maximal hotspot strength at the critical plasmon wavelength λ sp,c = L c . Figures 6(h) and (l) show maximal hotspot strength at λ sp,c ≈ 225 nm, which is in excellent agreement with the prediction of the geometric model. Furthermore, using m = 2, F 1 = 0.14 and, ρ = 0.23 (determined by exper imental measurements in the supplementary material) in  figure 4(a). The geometries (pink outlines) of the 1LG/2LG strip boundaries are revealed by (a) KPFM, while (b)-(f) s-SNOM amplitude images show the highly resonant nature of the plasmon interference hotspots that occur in the 1LG in that region. All scale bars are 1 µm. The white arrows in (a) indicate the scans of the (g)-(j) nearfield amplitude, s 2 , line profiles across the four numbered hotspot regions. The line profiles are normalised with respect to either the maximum recorded amplitude, s 2,max , or the graphene background, s 2,bg . Normalisation to the peak maximum in (g) is used to more clearly show the lateral movement of the hotspot (k)-(n) The spectral dependence of the hotspot locations and strengths, as extracted from the line profiles in (g)-(j), respectively. (k) The hotspot location, y , is defined in Figure 5(d). The hotspot gap in (m) is the spatial separation of the two hotspot peaks in (i). The dashed black lines are intended only as guides to the eye. equation (3b) gives maximum hotspot strength F HS,2 = F HS, + F 1 = 0.44 (i.e. hotspot strength of 1D cavity plus the first fringe from a single extra edge) for region 2, which is in good agreement with the maximum measured value of 0.43 taken from figure 6(l).
The following two regions, 3 and 4, have more complicated geometries and serve to test the limits of the model. Here, we approximate regions 3 and 4 as an equilateral triangle and a circle, respectively (see figure 5(c)). From equation (4), we predict λ sp,c = L c,3 / √ 3 for the triangle. Furthermore, we take λ sp,c = L c,∞ , where L c,∞ is defined as the critical diameter of the circle, as discussed in the previous section. As detailed in the supplementary material, we use KPFM measurements to give L c,3 ≈ 450 nm (λ sp,c ≈ 260 nm) and L c,∞ ≈ 240 nm (λ sp,c ≈ 240 nm) for the regions 3 and 4, respectively. A triangular constriction will behave as three wedges for λ sp < λ sp,c , which explains the behaviour seen in figures 6(b)-(f), where he hotspot in region 3 splits into three distinct hotspots as the incident wavelength λ 0 decreases. By taking a line profile across two of these hotspots (see figure 6(i)), we plot their separation as a function of plasmon wavelength, λ sp , (see figure 6(m)). Applying a linear fit, we find the critical plasmon wavelength, i.e. at which there is no gap between the hotspots, to be λ sp,c = 310 ± 60 nm, which is in reasonable agreement with the prediction of 240 nm. Finally, figures 6(j) and (n) show λ sp,c ≈ 210 nm for the hotspot in region 4, which is also consistent with geometric model of a subwavelength circle detailed above.
These results show that even approximating complex structures with simple geometric elements provides good quantitative predictions of the conditions of hotspot formation in sub-wavelength constrictions, as well as their strength. Higher resolution KPFM and AFM measurements would allow for the precision of this method to be pushed further.

Conclusions
Through the combination of high spectral resolution graphene plasmon nanoimaging and KPFM, we have validated a new wave model of surface polaritons in 2D nanoresonators. This model provides an intuitive understanding of complex phenomena using combinations of simple geometric approximations to provide quantitatively accurate predictions of hotspot behaviour. With this discovery, sophisticated graphene nanoresonators could be designed to manipulate the spectral and spatial nature of these hotspots, for use in plasmon-enhanced chemical analysis, environmental monitoring, and plasmonic nanoantennas for boosting the sensitivity of fluorescence microscopy and vibrational spectroscopy.