Introduction

Plasmons are optically driven collective electron oscillations capable of concentrating light down to deeply subwavelength volumes1, 2, laying the foundation for a multitude of essential applications in fields such as sensing3, 4, light-harvesting5, 6 and photochemistry7. While commonly used noble metals such as silver and gold exhibit plasmonic properties mostly in the visible and near-infrared region of the spectrum, graphene has recently emerged as a new material able to extend the field of plasmonics to infrared and terahertz wavelengths8, 9, 10, 11, 12, 13. Localized surface plasmon resonances can be excited in graphene by patterning it into optically resonant nanostructures, such as nanoribbons14, 15, 16, 17, 18, 19, 20 or nanodisks21, 22, 23. In contrast to noble metals, graphene plasmonic resonances can be dynamically tuned through electrostatic biasing24, 25 and enable a new generation of reconfigurable plasmonic devices26, 27, 28, 29, 30. In addition, plasmons excited in graphene can reach remarkably long lifetimes and provide unprecedented levels of light confinement31. Owing to these unique electro-optical properties, there has been a strong focus on developing new graphene-based plasmonic devices, and in particular infrared biosensors32, modulators33, 34, 35, metasurfaces36, 37, 38 and perfect absorbers39, 40.

In single-layer graphene (SLG), plasmonic resonances possess relatively low oscillator strengths because the attainable carrier densities are low. In the pursuit of stronger resonances, improved carrier mobility and new capabilities, the exploration of graphene plasmons in stack-ups composed of graphene and other two-dimensional (2D) materials has recently attracted significant interest. Initially, multilayer graphene (MLG) was investigated in the terahertz range and was shown to provide plasmonic resonances with larger amplitude modulation15. Extending this concept to infrared frequencies, more recent work has focused on bilayer graphene41, 42 and heterostructures composed of graphene on boron nitride31, 43. In addition, MLG structures have been shown to exhibit resonances with higher infrared transmission modulation compared to SLG44, 45. However, a satisfactory understanding of MLG, especially with regard to electrostatic biasing, is still lacking, thus limiting its application to actual devices. In this paper, we present a complete framework for understanding localized surface plasmon resonances in MLG stacks. In this multilayer structure, there is negligible interlayer electronic coupling and, unlike in graphite, the graphene band structure is retained. We explore the properties of MLG plasmonic nanoribbons, including spectral, near-field and tuning properties, and we show larger plasmonic enhancement associated with the multilayer geometry compared with previous results44. In particular, we show that double-layer graphene (DLG) provides wider spectral tuning ranges than SLG, in addition to sustaining resonances at higher frequencies and with higher extinction. Our analysis reveals that plasmon resonances in MLG can be quantitatively understood in the context of plasmon hybridization and provides a theoretical basis for the experimental evidence that the plasmonic response of DLG is equivalent to that of a SLG at a higher doping. Finally, we focus on the response of MLG to electrostatic biasing that demonstrate for the first time that the plasmonic effects of biasing are enhanced in MLG compared with SLG. Our results and proposed technique hold promise for a wide range of applications, as the performance of graphene-based photonic devices can be significantly improved by transferring two or more graphene layers rather than one.

Materials and methods

Fabrication and characterization

Graphene nanoribbon arrays are fabricated on a double-sided polished float-zone silicon substrate. A native oxide layer of 280 nm is grown by dry oxidation. Chemical vapor deposition (CVD)-grown graphene is wet-transferred over the oxide layer. DLG devices are formed by two successive transfers of SLG. Graphene nanoribbons are patterned with 100-keV electron beam lithography and are etched in oxygen plasma. The transmission spectra of the graphene arrays are acquired with a Fourier-transform infrared spectrometer coupled to an infrared microscope (Bruker Vertex 80 V and Hyperion 3000). An infrared polarizer is used to select the polarization transversal to the graphene nanoribbons. Measurements are performed in a dry-air-purged atmosphere with a non-immersion Cassegrain objective and a mercury–cadmium–tellurium detector cooled in liquid nitrogen.

Numerical simulations

The extinction spectra and charge distributions of the graphene nanoribbon arrays are calculated with a frequency-domain solver of Maxwell equations based on the finite elements method (commercial implementation Ansys HFSS). The periodic ribbon structure is simulated using Floquet–Bloch periodic boundary conditions. The convergence of simulations is ensured by monitorizing the transmission coefficient while refining the meshing in an iterative process. Graphene is modeled as a 2D layer with complex surface conductivity calculated from the Drude model: , where the Fermi level EF accounts for room temperature effects: . The frequency-dispersive dielectric function of SiO2 is obtained from Palik46.

Results and discussion

The excitation of localized surface plasmon resonances in a DLG nanoribbon array (2L-GNRA) is conceptually shown in Figure 1a. The nanoribbons have width W and period P. Scanning electron and atomic force microscope images of a typical 2L-GNRA (Figure 1b and 1c) show that the fabricated nanoribbons possess uniform width and well-defined edges. The infrared spectra of the 2L-GNRAs are presented in Figure 1d and are compared against those of single-layer nanoribbons (1L-GNRAs) with identical lateral dimensions. We observe that the extinction peaks produced by DLG have much higher intensity than those produced by SLG. The extinction intensity of 2L-GNRAs is approximately three times higher than that of 1L-GNRAs. This is a crucial advantage of DLG, which makes it very useful for infrared plasmonic applications where one of the main challenges is the poor coupling of free-space infrared radiation to the highly subwavelength and strongly confined plasmonic modes of graphene. Interestingly, unlike other techniques that rely on external resonant structures to increase graphene absorption47, the enhancement in DLG is intrinsic to the graphene layers and does not limit its frequency range of operation.

Figure 1
figure 1

Plasmonic resonances in SLG and DLG nanoribbons. (a) Conceptual representation of the excitation of localized surface plasmon resonances in a 2L-GNRA by illumination with mid-infrared radiation. The top graphene layer is transferred directly over the bottom layer without an interlayer spacer. Graphene nanoribbons are patterned with a width W and periodicity P. (b) Surface topography image measured with an atomic force microscope and (c) scanning electron microscope image of one of the fabricated 2L-GNRAs (W=50 nm, P=105 nm). (d) Experimental extinction spectra of 2L-GNRAs (solid curves) with different widths (W=25, 33, 40, 50 nm and P≈2W). The 2L-GNRA spectra are compared against those of 1L-GNRAs (dashed curves) with identical width and periodicity. The extinction is calculated as 1−T/T0, where T and T0 are the transmission coefficients with and without GNRAs, respectively.

Full size image

In addition to the higher extinction, we observe that 2L-GNRAs resonate at higher frequencies with respect to 1L-GNRAs. Therefore, it is possible to obtain a plasmonic resonance at a specific frequency (wavelength) using wider DLG ribbons instead of narrower SLG ribbons. This is particularly useful from a fabrication perspective because the nanostructuring of ribbons with highly uniform width poses significant challenges, especially at infrared frequencies for which the required ribbon widths are of the order of tens of nanometers. A further important advantage is that the quality factor of the plasmonic resonances in DLG is essentially maintained with respect to SLG. This indicates that the carrier scattering rate and mobility in DLG and SLG are similar.

The spectra of 2L-GNRAs do not present any sharp feature superimposed on the plasmonic resonance, indicating that DLG is free of infrared-active phonons in this spectral range. For most applications, such as plasmon-enhanced biosensing, these additional features in the plasmonic response are undesired32. In this sense, DLG presents a clear advantage over other related graphene materials such as bilayer graphene or graphene–boron nitride heterostructures, which exhibit infrared-active phonons at 1580 and 1370 cm−1, respectively41, 43, 48. We remark that the DLG used in this work should not be confused with bilayer graphene, in which the atomic alignment of the two graphene layers strongly modifies the electronic band structure and ensues optical properties41, 49. Because of our fabrication approach, in which we perform two independent transfers of CVD graphene, the interlayer separation is larger than in bilayer graphene and there is no crystallographic alignment; therefore, the plasmonic effects of interlayer tunneling are negligible42. Roughness measurements of the bottom graphene layer indicate a minimum average interlayer spacing of 1.0 nm because of the presence of polymeric residues and graphene wrinkles.

We study next the optical response of the GNRAs when the Fermi level of graphene is electrostatically tuned. The tuning scheme is outlined in Figure 2a and is based on applying a voltage difference across the SiO2 layer to inject or remove electrons from the graphene layers and modify their Fermi level. In contrast to previous work15, our stack geometry avoids using isolating spacers between layers and keeps the graphene layers electrically interconnected, allowing them to be simultaneously biased. The physical mechanisms underlying the biasing of MLG are analyzed later in this manuscript. The voltage-dependent spectra are shown in Figure 2b and 2c for 1L- and 2L-GNRAs, respectively. For 1L-GNRAs, we observe that biasing from 0 V to negative voltages (up to −100 V) produces a progressive shift of the plasmonic resonance toward higher frequencies and an increase in the resonance intensity. This result is consistent with previous experiments reported by other research groups14, 16, 21 and is a consequence of the initial p-doping of our CVD graphene (EF0<0), for which negative voltages produce an increase in the hole carrier density. The results of biasing 2L-GNRAs with negative biasing voltages show identical behavior as in 1L-GNRAs, that is, a progressive blue shift of the resonance accompanied by an increase in its extinction. A clear advantage of 2L-GNRAs is that for the maximum negative bias they reach higher frequencies and intensity than 1L-GNRAs with the same lateral dimensions. In addition, while in 1L-GNRAs the use of positive bias voltage quenches the plasmonic resonance, this does not happen in 2L-GNRAs where the much higher extinction of the plasmonic peaks makes it possible to use positive biasing. Figure 2c shows that positive bias voltages push the plasmonic resonance toward lower frequencies and thus increase the spectral tuning range.

Figure 2
figure 2

Electrostatic tuning of plasmonic resonances in 1L- and 2L-GNRAs. (a) Electrostatic biasing scheme of graphene nanoribbon arrays. The biasing voltage is applied across the 280-nm-thick SiO2 layer. Negative (positive) voltages shift the Fermi level (EF) of our p-doped graphene away from (toward) the Dirac point. (b, c) Experimental extinction spectra of 1L- and 2L-GNRAs for different biasing voltages Vg. For 1L-GNRAs the voltage Vg is varied between 0 V (black) and −100 V (red), whereas for 2L-GNRAs Vg is varied between +100 V (blue) and −100 V (red). (d, e) Extinction value and spectral position of the measured plasmonic peak in 1L- and 2L-GNRAs for different ribbon widths W and biasing voltages Vg.

Full size image

We perform an extended comparison between SLG and DLG in Figure 2d and 2e, where we present the spectral position and extinction of the plasmonic peaks for all the fabricated arrays and tested bias voltages. Each marker in the plot corresponds to a specific array and voltage (Vg), with its horizontal position indicating the resonance frequency and the vertical axis showing the intensity of the resonance peak. At least two identical arrays were fabricated for each value of the width W, which correspond to the overlapping lines in the graph, demonstrating excellent repeatability. For a variable voltage Vg there is a linear relationship between the resonance peak intensity and spectral position for both 1L- and 2L-GNRAs. The linear relationship is only broken at the lowest frequencies, where SiO2 is strongly dispersive and its effect is equivalent to a compression of the frequency axis. The comparison between Figure 2d and 2e highlights very clearly the advantages of DLG over SLG: higher peak intensity (up to 10% for 2L and 6% for 1L), higher resonance frequency (up to 1950 cm−1 for 2L and 1800 cm−1 for 1L) and wider tuning range (up to 450 cm−1 for 2L and 300 cm−1 for 1L).

To understand the physical mechanisms controlling the plasmonic response of DLG, we perform numerical simulations for a 2L-GNRA composed of two graphene layers separated by a distance s and with Fermi levels EF(1) and EF(2) defined relative to their respective Dirac points. In the spectral range of interest, the conductivity of graphene follows the Drude model50, 51

Figure 3a shows the evolution of the simulated 2L-GNRA extinction spectra with interlayer separation s. In these calculations, we study plasmonic effects for non-varying graphene parameters and exclude, for instance, potential variations of Fermi level induced by changes in the interlayer distance. For large separations (s>W) the 2L-GNRA produces two resonances matching those of the individual layers, that is, the resonance frequency of 1L-GNRAs with Fermi levels EF(1) and EF(2) (dashed lines). This indicates that for a separation comparable to or larger than the ribbon width, the two graphene layers are electromagnetically uncoupled. As we progressively reduce the interlayer separation s, there is an increasing spatial overlap between the near-fields of the two graphene layers and therefore higher electromagnetic coupling52. This is evidenced in our simulations by the progressive spectral shift of the two resonances as the separation s approaches zero. The simulated plasmonic response is in perfect agreement with a hybridization model that predicts symmetric and antisymmetric modes with resonance frequencies above and below those from the two individual layers, respectively53, 54, 55. The simulated charge distribution corresponding to these two modes is shown in Figure 3b. It is clear that the dipolar modes of each layer have the same polarity for the symmetric mode and opposite polarity for the antisymmetric mode. When the separation s becomes much smaller than the width W, the antisymmetric mode is extinguished (that is, it becomes dark) and only the symmetric mode remains. It is this symmetric mode that is detected in our 2L-GNRAs experimental results, for which we measure a higher resonance frequency and extinction than for 1L-GNRAs.

Figure 3
figure 3

Physical mechanisms and hybridization in 2L-GNRAs. (a) Simulated extinction coefficient of a 2L-GNRA for a variable interlayer separation s. The two graphene layers have Fermi levels |EF(1)|=0.3 eV and |EF(2)|=0.2 eV, relaxation time τ=100 fs, the width of ribbons is W=40 nm and are embedded in a dielectric with refractive index n=1.4. Solid curves represent the resonance frequencies of the 2L-GNRA and dashed curves stand for the resonance of each individual layer, which correspond to the resonance of a 1L-GNRA with Fermi level |EF(1)| or |EF(2)|. (b) Electric charge distribution across the ribbons for the 2L-GNRA with s=W/4 at the two resonance frequencies indicated in the figure. The high- (low-) resonance corresponds to a symmetric (antisymmetric) mode. (c) Simulated extinction spectra of the 2L-GNRA (solid curves) for an interlayer separation varying from s=W (red curve) to s=W/100 (blue curve). Dashed curves represent the extinction spectra of 1L-GNRAs with Fermi levels |EF(1)|, |EF(2)| and |EF(1)|+|EF(2)|. (d) Field enhancement distribution at resonance frequency for 1L-GNRA (|EF|=0.3 eV) and for the symmetric mode of a 2L-GNRA (|EF(1)|=0.3 eV, |EF(2)|=0.2 eV) with a separation s=1 nm. The field is presented for the cross-section indicated in the inset. (e) Field enhancement for the 1L- and 2L-GNRAs along the white dashed line in subfigure d.

Full size image

The response of the 2L-GNRA is further analyzed in Figure 3c, where we compare the extinction spectra of 1L-GNRAs (dashed lines) and 2L-GNRAs (solid lines) for an interlayer separation s varying from s=W (red curve) to s=W/100 (blue curve). As expected, when the two layers are uncoupled (s=W) there is nearly perfect overlap between the 2L-GNRA spectrum and that of two independent 1L-GNRAs with Fermi levels EF(1) and EF(2). Interestingly, when the two layers are fully coupled (s<<W) the 2L-GNRA spectrum perfectly matches the spectrum of a 1L-GNRA with an equivalent Fermi level EF2L, given by the sum of the Fermi levels of the individual layers (EF2L=|EF(1)|+|EF(2)|)56, 57, 58, 59. At first glance, this result is counterintuitive in the context of collective electronic excitations in coupled 2D systems60, since one might expect the equivalent carrier density of the two-layer system to be the sum of the carrier densities in each layer (ns(1)+ns(2)), which is not the case here. The explanation for the equivalence between a 2L-GNRA with Fermi levels EF(1), EF(2) and a 1L-GNRA with Fermi level |EF(1)|+|EF(2)| is that for small separations one can replace the two parallel conductive layers by an equivalent layer having the sum of the two conductivities. For the Drude-model graphene conductivity Equation (1)), the resulting equivalent conductivity of the two-layer system is

This expression corresponds to the Drude conductivity for a graphene layer with Fermi level

Accordingly, we define the effective doping of DLG as the sum of the individual doping levels of its two constituent layers. Taking into account the relation between Fermi level and carrier density (), the equivalent carrier density of DLG is

It is important to note that, even though the present work focuses on the nanoribbon geometry, the results are equally valid for any arbitrary shape, as long as the separation between the layers is much smaller than the plasmon wavelength and the graphene conductivity follows a Drude model. In this model, the plasmon frequency scales as with the Fermi energy, and, consequently, the 2L-GNRA has a higher resonance frequency than a 1L-GNRA with the same total carrier density (ns(1)+ns(2)), reaching a frequency increase factor of 21/4 (that is, 20% higher) when the carriers are equally distributed among the two layers.

This hybridization model can easily be extended to a higher number of layers. An analogous study for a three-layer graphene structure (Supplementary Fig. S1) shows that the different polarities of the dipolar modes on each graphene layer give rise to multiple modes. As in the 2L-GNRA, when the interlayer separation s approaches zero, all of the modes become dark, except the so-called bonding mode in which all layers have the same polarity. For small separations (s≪W) the spectrum of the multilayer structure matches that of a 1L-GNRA, adding up the Fermi levels of all layers () or equivalently having a carrier density . Therefore, increasing the number of layers is a practical approach to create graphene structures with higher equivalent doping. However, a word of caution has to be given: the higher effective doping of DLG does not extend the spectral region that is free from interband transitions ()50, a condition that is required in order to have long-lived plasmons that are unaffected by Landau damping.

A study of the plasmonic field distributions is presented in Figure 3d, where we compare the near fields in 1L- and 2L-GNRAs. We observe that both structures display field distributions corresponding to dipolar modes with intense near-field hotspots at the graphene nanoribbon edges. Furthermore, in Figure 3e we present the field intensity profile perpendicularly to the nanoribbons close to the ribbon edge (0.3 nm distance), thus avoiding any inaccuracy introduced by classical electromagnetic computations near sharp boundaries. We observe that the field distribution, confinement and optical enhancement in DLG do not differ significantly from those in SLG.

The plasmonic equivalence between SLG and DLG can be experimentally demonstrated by comparing the measured spectra of the respective nanoribbon arrays. Figure 4a shows the extinction spectra of 1L-GNRA and 2L-GNRA for specific pairs of biasing voltages. We can clearly see that the spectrum of the 2L-GNRA for a given voltage is in excellent agreement with that of the 1L-GNRA for a significantly lower (more negative) voltage. These results provide experimental confirmation that a 2L-GNRA is equivalent to a 1L-GNRA with a different Fermi level. The change of the equivalent Fermi level with bias voltage in MLG is investigated next.

Figure 4
figure 4

Equivalent doping and enhanced tunability of DLG. (a) Comparison between the experimental extinction spectra of 1L- and 2L-GNRAs with W=33 nm and specific pairs of bias voltage Vg. (b) Band diagram representation of a multilayer graphene stack composed of two layers that are electrically interconnected. The black/red planes represent the Fermi level before/after biasing. (c) Band diagram representation for initially undoped graphene layers being biased. The Fermi levels and carrier densities are represented for SLG (1L), DLG (2L) and its two individual layers (1) and (2). (d) Simulated extinction spectra of the 1L-GNRA with W=33 nm, graphene relaxation time τ=15 fs and varying Fermi level EF. (e) Equivalent carrier density and Fermi level of the 1L- and 2L-GNRA extracted from experimental spectra. The dashed and solid curves represent a linear fit of the equivalent carrier density for the 1L- and 2L-GNRA, respectively.

Full size image

We consider a stack of N graphene layers having different levels of initial doping ( for the i-th layer), being electrically interconnected to be simultaneously biased and having an interlayer separation s tending to zero. In this case, the Dirac cones of the interconnected layers shift along the energy axis; therefore, the Fermi levels of all layers are aligned61, 62, 63, 64, as shown in Figure 4b for two layers. Increasing the bias voltage Vg shifts the overall Fermi level and as a result the injected carriers are redistributed among the different layers proportionally to their doping level. The Fermi level of each layer after applying a bias () is calculated by assuming that the Fermi-level variation is the same for all layers and that the extra carriers in the graphene layers () add up to the total charge in a capacitor with surface capacitance Cox:

where the Fermi level of each layer is defined as relative to its corresponding Dirac point. This set of second-order equations () has analytical solution in the case where all graphene layers have the same type of doping (p or n). After some calculations, the variation of Fermi level (ΔEF in Figure 4b) can be written in the simple form (see the Supplementary Information for a detailed derivation)

where and the plus/minus sign accounts for the doping type of the graphene layers (p- or n-doped). Notice that the Fermi level variation does not depend on the individual doping levels of the graphene layers (), as one might initially expect. Instead, it depends only on the sum of the initial levels (), which we have previously shown to correspond to the equivalent Fermi level of the MLG. Strictly speaking, the equivalent Fermi level corresponds to the sum of the absolute doping levels (); however, the expression without absolute values can be used when all layers share a common doping type. Next, we calculate the equivalent Fermi level of the biased MLG () by adding Equation (7) for all layers

or in terms of the equivalent carrier density ,

This simple result shows that the increase in the equivalent carrier density of MLG is enhanced with respect to SLG by a factor corresponding to the number of layers N. This insight is of crucial importance for plasmonic applications, since it shows that MLG provides broader tunability than SLG. In particular, the biasing of DLG can be twice as effective as in SLG, meaning a twofold increase in the equivalent carrier density for the same bias voltage. This effect can be readily understood for the case of undoped graphene layers (), where the carriers split equally between the two layers () and from Equation (4) the equivalent DLG doping is (Figure 4c). Our simple model provides a clear and intuitive description of the effects of electrostatic tuning over the plasmonic response of MLG; however, more elaborate theoretical models may be required for accurate calculations of charge redistribution in structures with a large number of graphene layers65, 66, 67. The biasing effectiveness can be degraded if carriers are redistributed differently between the two graphene layers. However, the impact of a moderate charge redistribution over the equivalent carrier density is relatively low. For instance, if an undoped layer receives four times as many carriers as a second undoped layer, the equivalent carrier density according to Equation (4) is reduced by just 10% with respect to an even distribution of charges.

To experimentally validate the enhanced tuning capabilities of MLG, we extract the equivalent Fermi level of our 2L-GNRA for the different voltages and compare it against that of the 1L-GNRA. Specifically, we extract the Fermi level by comparison with electromagnetic simulations. In Figure 4d we show the simulated extinction spectra of our 1L-GNRAs for different Fermi levels. We observe that the spectral lineshapes are in good agreement with the experimental data. Since the resonance frequency shifts monotonously with the Fermi level, we can use the simulations to connect these two parameters and extract the Fermi level from the experimental data. Figure 4e presents the extracted carrier densities and Fermi levels of SLG (EF1L, ns1L) and DLG (, ). As expected, the equivalent doping level in DLG is significantly higher than in SLG, confirming that DLG is equivalent to highly doped SLG. We observe two distinct behaviors in the DLG carrier density for voltages Vg above or below +40 V. For high positive voltages, we observe a saturation of the DLG equivalent carrier density. This saturation is caused by the Fermi level (red plane in Figure 4b) crossing the Dirac point of one of the layers, resulting in two layers with opposite doping types (see Supplementary Fig. S3).

Next, we focus on voltages below +40 V, which correspond to both graphene layers being p-doped. This is the relevant voltage range for plasmonics, since it corresponds to a high doping. In this range, the equivalent carrier density varies linearly with the applied voltage, in agreement with our theoretical calculations (see Equation (7)). This linear relationship is well known for SLG, but is observed here for the first time in DLG. More importantly, the variation of the equivalent carrier density with the bias voltage (dns/dVg) in DLG is 36% higher than in SLG. This variation, despite being smaller than the theoretical value, confirms that biasing is more effective in MLG structures than in SLG. It is important to emphasize that the improved biasing is not caused by a larger amount of total injected carriers (which depends only on the bias voltage and the gating layer capacitance); instead, it is due to the redistribution of these carriers over multiple plasmonically coupled layers. Finally, it is clear that stacking multiple graphene layers is a powerful approach to enhance the tunability of graphene plasmonic devices.

Conclusions

In conclusion, we have shown that DLG, readily fabricated by two successive transfers of CVD graphene, provides significant advantages over SLG for plasmonic applications. DLG operates as a SLG with an equivalent higher doping and serves as a plasmonic platform for higher frequency resonances, with stronger intensity and wider tuning ranges. These results can be further generalized to MLG structures with more than two layers, potentially yielding even better performances. Since high doping and broad tunability are critical requirements for graphene to succeed in plasmonic applications, we believe that the combination of multiple graphene layers opens very promising pathways toward novel reconfigurable nanophotonic devices.