Tunable plasmons in regular planar arrays of graphene nanoribbons with armchair and zigzag-shaped edges

Recent experimental evidence for and the theoretical confirmation of tunable edge plasmons and surface plasmons in graphene nanoribbons have opened up new opportunities to scrutinize the main geometric and conformation factors, which can be used to modulate these collective modes in the infrared-to-terahertz frequency band. Here, we show how the extrinsic plasmon structure of regular planar arrays of graphene nanoribbons, with perfectly symmetric edges, is influenced by the width, chirality and unit-cell length of each ribbon, as well as the in-plane vacuum distance between two contiguous ribbons. Our predictions, based on time-dependent density functional theory, in the random phase approximation, are expected to be of immediate help for measurements of plasmonic features in nanoscale architectures of nanoribbon devices.


Introduction
Quantized, coherent and collective density fluctuations of the valence electrons in low-dimensional nanostructures, better known as plasmons, have been attracting significant interest, due their capability to couple with light and other charged particles, thus paving the way to novel applications in a wide range of technologies, as diverse as biosensing, light harvesting or quantum information [1][2][3][4][5]. On more fundamental grounds, plasmon-like modes are the "true" low-energy excitations of low-dimensional systems [6,7], while charged and spinful modes are realized as coherent states, with their own peculiar dynamics [8,9], both in normal superconducting phases [10][11][12][13][14]. Graphene has first emerged as an extraordinary platform for controlling the propagation of surface-plasmon waves [15], because of its unique electronic and optical properties [16]. In particular, the extrinsic plasmons of this one-atom-thick hexagonal lattice of sp-bonded carbon atoms have shown much stronger confinement, larger tunablity and lower losses with respect to more conventional plasmonic nanoparticles, such as, for example, silver and gold [17]. With the rise of low-dimensional materials, a number of theoretical and experimental studies have been oriented to launch, control, manipulate and detect plasmons in graphene-related and beyond-graphene structures [18][19][20][21], which are expected to be embedded in nextgeneration nano-devices that may operate from infrared (IR) to terahertz (THz) frequencies [22][23][24][25]. As a noteworthy example, graphene nanoribbons (GNRs) preserve most of the exceptional features of graphene, with the additional property that they are semiconductors and their band gap is geometrically controllable.
A clear picture of confined edge (interband) and surface (intraband) plasmons in GNRs, as wide as 100-500 nm, has been achieved by infrared imaging measurements on the nanoscale [26]. On the theoretical side, some nearest-neighbor tightbinding [27,28] and semiclassical electromagnetic [29,30] approaches have been able to characterize the intraband mode, being generally excited by a THz electromagnetic field pulse. Very recently, an ab initio analysis has elucidated the role of both intraband and interband plasmons in narrow GNRs below ca. 1 nm in width [31].
In this work, we extend the results of [31], to include some GNRs up to 2.2 nm wide, and demonstrate how the extrinsic plasmons of these systems can be finely tuned by changing a small number of simple parameters, such as the unit-cell length, width and chirality of each GNR, or the in-plane distance between contiguous GNRs. Our study is carried out using time-dependent (TD) density functional theory (DFT), within the random phase approximation (RPA). The computations are performed at room temperature (T = 300 K), including both intrinsic and extrinsic conditions. Different zigzag (Z) and armchair (A) configurations are examined, with the GNR-ends being passivated by hydrogen atoms, which mimics an ideal setup of long ribbons, with perfectly symmetric edges, suspended or grown on inert substrates. A specific focus is made on the 5AGNR, 11AGNR and 4ZGNR, 10ZGNR geometries that are, respectively, characterized by 4,10 zigzag-chains and 5,11 dimer-lines across the width of the GNRs. The dielectric properties of these systems are calculated in response to probe electrons or photons with incident energies, ω, smaller than 20 eV, and in-plane incident momentum q below 0.8 Å −1 . For comparison purposes, the well-known intrinsic plasmonics of graphene are also reported.
In the following, we briefly account for the theoretical tools that we have used to explore the electronic structure and dielectric properties of 5AGNR, 11AGNR and 4ZGNR, 10ZGNR and graphene (Results and Discussion, chapter 'Theoretical framework'). Next, we will present a detailed analysis of the changes induced by the above mentioned parameters to the peak position, intensity, and dispersion of the GNR-plasmons (Results and Discussion, chapter 'Tunable plasmons in GNR arrays').

Theoretical framework
Our TDDFT approach is divided into two steps. First, the (ground-state) electronic properties of the different GNRs (and graphene) are obtained by DFT. Second, the basic equation of linear-response theory in the RPA is employed, with a corrected electron-electron interaction, to calculate the dielectric properties, and hence the plasmon structure, of the systems.

DFT method
Density-functional calculations are performed using the planewave (PW) basis-set [32], i.e., the normalized space functions , which depend on the wavevectors k of the first Brillouin zone (1st BZ), and the reciprocal lattice vectors G associated to the three-dimensional (3D) crystal of unit-cell volume Ω 0 . The ground-state electronic properties of the different GNRs (as well as graphene) are computed within the local density approximation (LDA), being defined by the Perdew-Zunger parameterization of the uniformgas correlation energy [33]. Norm-conserving pseudopotentials of the Troullier-Martins type are adopted to eliminate the core electrons [34]. A cut-off energy of ca. 680 eV on the number of PWs is sufficient to obtain well-converged electronic structures with ca. 10 5 PWs.
In the GNR arrays, the C-C length is allowed to range from 1.414 to 1.426 Å, while the C-H bond length is fixed to 1.09 Å, with a bond-angle of 120 °(as shown in [31] relaxation effects play a minor role in 4ZGNR and 5AGNR). The 3D periodicity required by PW-DFT is generated by using an in-plane vacuum distance of 5-20 Å, and an out-of-plane lattice constant of 15 Å. The self-consistent run is carried out using an unshifted (and Γ-centered) Monkhorst-Pack (MP) grid, made of 60 × 60 × 1 k-points [35], which results in a uniform sampling of the irreducible 1st BZ on the ΓX-segment (shown in Figure 1). The converged electron densities are subsequently used to compute the Kohn-Sham (KS) electronic structure on a denser MP-mesh of 180 × 1 × 1 k-points, including up to 120 bands, which is enough to explore the dielectric properties of the GNRs at probing energies below ≈20 eV. The IR to THz region is further scrutinized with a finer MP mesh of 2000 × 1 × 1 k-points, including up to 30 bands. The main results of our DFT computations are summarized in the plots of Figure 1, which show the different geometry, band structure and density of states (DOS) of the GNR arrays. 4ZGNR and 10ZGNR behave as semimetals, with barely touching valence and conduction bands (Figure 1a,c). 5AGNR and 11AGNR are small-gap semiconductors (Figure 1b,d), contrary to nearest-neighbor TB approaches in which all AGNRs appear gapless [27].
Indeed, several DFT studies have carefully characterized the band gaps of ZGNRs and AGNRs [36][37][38][39]. In particular, local spin-density calculations suggest the opening of a band gap larger than 0.1 eV in ZGNRs [36,37]. However, the LDA analysis of a virtually gapless GNR, i.e., 4ZGNR or 10ZGNR, in comparison with 5AGNR and 11AGNR, is particularly instructive to emphasize the different role played by doping in separating the extrinsic plasmon modes, which will be detailed in chapter 'Tunable plasmons in GNR arrays'.
The simulations at hand explore ranges of different geometrical and conformation parameters, which will be used to characterize the tunability of the GNR plasmons. In particular, (i) GNR widths (w) of around 0.7-2.2 nm are sorted out; (ii) zigzag and armchair edges are considered, to elucidate the role played by chirality; (iii) in-plane vacuum distances from 5 to 20 Å are tested; and (iv) different unit-cell extensions are simulated by changing the C-C bond length by about 0.5%, to account for stretching effects. As for intrinsic graphene, the C-C bond length and out-of-plane lattice constant are fixed to 1.42 Å and 15 Å, respectively. The self-consistent run is performed on a 60 × 60 × 1 MP-grid, and the KS electronic properties are then computed on a 180 × 180 × 1 MP-grid, including up to 80 bands.

TDDFT approach
The KS eigenvalues (see Figure 1) and corresponding eigenfunctions are the main outputs of the DFT computations, with N being the total number of k-points in the 1st BZ, and ν being the band index. The KS eigensystem gives access to the unperturbed density-density response function, of the non-interacting valence electrons, to a probe particle of energy ω and momentum q. The latter is provided by the Alder-Wiser formula [40,41]: (1) Indeed, plasmons in solid-state materials are typically triggered by electron-beam radiation, photo-currents, and even charged ions, with incident kinetic energies of the order of 0.1-1 keV [42,43]. In the present context, the probe particle is an electron or a photon that weakly perturbs the system.
In Equation 1 (expressed in Hartree atomic units), the factor of two takes into account the electron spin, η indicates an infinitesimal (positive) broadening parameter (set to 0.02 eV), f νk is the Fermi-Dirac distribution, and (2) labels the density-density correlation matrix elements, which depend on the number of PWs included in the DFT simulations.
The full susceptibility or interacting density-density response function is determined by the central equation of linearresponse TDDFT [44,45] (3) where, in the RPA, the υ GG′ terms are approximated to the bare Coulomb coefficients: (4) However, the long-range character of the Coulomb potential yields non-negligible interactions between replicas along outof-plane direction z. To erase this unphysical phenomenon, a two-dimensional (2D) cut-off, say, a truncated Coulomb potential is used [46][47][48][49] to replace the υ GG′ terms by the truncated Fourier integral: (5) Here, g and G z denote the in-plane and out-of-plane components of G, respectively.
With this 2D correction in mind, we can introduce the inverse dielectric matrix: (6) The zeros in the real part of the macroscopic dielectric function (permittivity) provide the condition for a plasmon resonance to occur, stated as: (7) The imaginary part of the inverse permittivity is proportional to so-called energy loss (EL) function, which provides the plasmon structure: Non-local field effects [50] are included in E LOSS through Equation 3. We have verified that ca. 120 G vectors for all GNRs (and ca. 51 G vectors for graphene), sorted in length and being of the form (0,0,G z ), give well-resolved and converged results in the sampled energy-momentum region, delimited by ω < 20 eV and q < 0.8 Å −1 .

Tunable plasmons in GNR arrays
We proceed by clarifying the role played by the geometric and conformation parameters introduced above in the different GNRs, whose intrinsic response is shown in Figure 2 together with that of graphene. In the following, we will also evaluate several extrinsic conditions associated to Fermi energy shifts ΔE F in the range of −0.2 to 0.2 eV.

Ribbon width and chirality
Independently on the width and chirality of the ribbons, all GNRs are characterized by two interband plasmons at excitation energies above about 2 eV. These excitations, shown in Figure 2b-e, are analogous to the well-known π and π-σ plasmons of graphene [51,52], as displayed in Figure 2a. Similar features occur in bilayer graphene [51], multilayer graphene is represented as a density-color plot vs incident energy ω (in eV) and momentum q (in Å −1 ). The latter is sampled along the ΓΚ path of the 1st BZ of graphene (f), and the ΓΧ path of the 1st BZ of the different GNRs (g). Besides the π and π-σ plasmons, the ZGNRs present a low-energy intraband mode (IntraP), whereas the AGNRs have a low-energy interband mode (InterP). The intensity scale in (a-d) is cut at 80% of the π peak maximum, to ease comparison between the different density-color plots.
The energy window displayed in Figure 2 does not show the complete energy-momentum dispersion of the π-σ plasmon, however, the latter seems to be quadratic in graphene and linear in the GNRs. At long wavelengths, in the q→0 region, the π plasmon of all systems has a -like dispersion, while at q > 0.2 Å −1 it exhibits linear behavior. The intrinsic plasmons of 10ZGNR and 11AGNR appear in the same energy region as those of graphene, i.e., at ω ≈ 4-5 eV and ω ≈ 14-15 eV. On the other hand, they are red-shifted in 4ZGNR and 5AGNR, with the π plasmon having a peak at ω ≈ 2-3 eV and the π-σ plasmon lying at ω ≈ 13-14 eV. Furthermore, the π and π-σ plasmons detected in 4ZGNRs (w ≈ 0.9 nm) and 5AGNR (w ≈ 0.7 nm) exhibit markedly discontinuous dispersions, being split into more branches [31]. This is a consequence of the narrow widths of the two systems that generate several, distinct one-dimensional bands of π-and σ-character (Figure 1c,d).
By increasing the GNR width (w > 1 nm), the number of bands increases, and less disjoint plasmon dispersions appear, which clearly tend to the continuous patterns of graphene (w→∞). Thus, semiconducting and semimetallic GNRs have plasmon resonances in the visible (VIS) to ultraviolet (UV) regime that may be controlled by the GNR width. This tunability feature is evidently absent in graphene.
Quantum confinement and chirality are key factors for plasmon resonances at frequencies smaller than 2 eV. We see that zigzag systems exhibit an intraband plasmon, while armchair systems present an interband plasmon. These two modes correspond to the surface and edge plasmons detected in large-width, extrinsic GNR arrays fabricated on Al 2 O 3 [26]. The surface plasmon of ZGNRs is originated by the large DOS peak observed at the Fermi level E F (Figure 1a,c). This mode shows a -like dispersion [31] and seems to be analogous to the conventional 2D plasmons of extrinsic graphene [19,51]. The edge plasmon of AGNRs appears as an effect of collective excitations generated close to E F [31], associated to single-particle excitations that connect the two DOS peaks around E F (Figure 1b,d). The characteristics of this interband mode are similar to those of the π plasmon of intrinsic graphene, i.e., at long wavelengths the interband plasmon shows a -like dispersion, while at q > 0.1 Å −1 it displays a linear dispersion.
Both the intraband and interband modes have been proved to be genuine plasmons in intrinsic 4ZGNR and 5AGNR [31], as they are strictly associated to the zeroes of the real permittivity, satisfying the condition given by Equation 7. It has been further demonstrated that extrinsic 4ZGNR presents only an intraband plasmon structure, independently on the positive doping level used (below ca. 1 eV), while both intraband and interband plasmons coexist in 5AGNR [31].
To support this result, we report in Figure 3 the macroscopic dielectric function and the EL function of the different GNR arrays for a selected momentum value (q = 0.011 Å −1 ) and a negative doping level (ΔE F = −0.1 eV). We see that 10ZGNR and 4ZGNR present similar plasmonic features, with the intraband plasmon resonance being blue-shifted by increasing the GNR width (Figure 3a,c). In 11AGNR and 5AGNR, not only the peak position but also the interplay of the interband and intraband plasmon is strongly dictated by the doping level and the GNR width (Figure 3b,d). In 5AGNR, the two modes are well resolved in energy, with the zeroes of the real permittivity being hidden by the Landau damping mechanism, associated to single-particle excitation processes [25,[46][47][48]. In 11AGNR the same modes strongly interfere and largely dominate with respect to single-particle excitations. A similar interplay was observed in extrinsic 5AGNR subject to a positive doping of about 0.3 eV [31]. These outcomes are basically due to the different band-gap values of the two AGNRs, which according to our predictions are ca. 0.18 eV for 11AGNR, and ca. 0.36 eV for 5AGNR. Accordingly, less energy requirements are needed to produce a welldefined intraband collective electronic excitation in 11AGNR. On the other hand, a positive doping larger than 0.2 eV yields a well-defined intraband plasmon in 5AGNR [31]. Interestingly enough, some GNRs with band-gap values of the same order of 11AGNR and 5AGNR have been recently synthetized on Au(111) [60]. Then, our ab initio analysis can be of help in interpreting plasmon measurements on currently synthetized GNR-structures.
Chirality seems to be a major point for the design of GNRbased plasmonic devices. One or two plasmon modes can be exploited, depending on the shape of the GNR edges. In this respect, negative or positive doping acts as a modulating factor of the plasmon modes.
In Figure 4 we see that a change in doping sign, from −0.1 to 0.1 eV, produces a slight red shift in the intraband plasmon of 10ZGNR and the interband plasmon of 11AGNR (Figure 4a,b). More significant variations are observed in the intraband plasmon of 11AGNR, which is markedly blue-shifted and doubled in intensity by the same change of extrinsic conditions (Figure 4b). Therefore, an asymmetric response is observed in the intraband plasmon of semiconducting GNRs (Figure 4b). Moreover, as the GNR width decreases an appreciable blue/red shift is detected in the plasmon peaks of both ZGNRs and AGNRs (Figure 4c,d). Thus, a tunable energy response may be more strongly influenced by the ribbon width than the doping level. In (a-c) the energy-momentum region ω ≤ 1.6 eV and q || ΓX ≤ 0.1 Å −1 is explored, showing how the parameter L modifies the relative position between the intraband and interband plasmons, denoted IntraP and InterP as in Figures 2-4. The intensity scale is cut at 95% of the IntraP peak. The effect of changing the in-plane vacuum distance is even more evident in (d), where the different EL functions of (a-c) are compared at a fixed momentum value of q = 0.025 Å −1 parallel to ΓΧ.

Mechanical deformations
Let us now see how the fascinating plasmonic features of semiconducting GNRs are affected by changes of the in-plane separation. With reference to the case of 5AGNR, we take a positive doping value of 0.2 eV and consider vacuum distances L, between continuous arrays, in the range of 5 to 20 Å ( Figure 5).
As a first result, we see that both intraband and interband plasmon modes exist in 5AGNR, no matter how far apart the arrays are. The intraband plasmon is, however, affected in intensity, while the interband plasmon is blue-shifted as the vacuum distance decreases down to 5 Å. This effect is clearly visible at q = 0.025 Å −1 in Figure 5d, where a broad interband plasmon peak is detected at ω ≈ 0.6-1 eV. Both the large blue shift of the interband plasmon, and the intensity decrease of the intraband plasmon, are consistent with the idea that as the GNR arrays get closer a large graphene area is created.  In (a-c), the energy-momentum region ω ≤ 1 eV and q || ΓX ≤ 0.03 Å −1 is explored, showing that a cc (or the lattice constant 3a cc ) is a major factor in modulating the intraband and interband plasmons, denoted IntraP and InterP as in Figures 2-5. The intensity scale is cut at 40% of the IntraP peak. The effect of stretching the 5AGNR structure is even more evident in (d), where the different EL functions of (a-c) are compared at a fixed momentum value of q = 0.025 Å −1 parallel to ΓΧ.
When the vacuum distance becomes negligibly small, the interband plasmon detected in AGNRs enters the region where the π plasmon of graphene are found, while the intraband plasmon decreases in intensity to a small contribution, reported in room temperature calculations of slightly doped graphene [19,31].
Finally, we show how the intraband and interband plasmons of 5AGNR are affected by stretching/shrinking the unit cell of the system by about 0.5%, with respect to its nominal value associated to a C-C bond length, a, of 1.42 Å. In this application, the in-plane vacuum distance is fixed to 15 Å is and a negative doping level of −0.2 eV is considered. As shown in Figure 6, the band gap decreases with increasingly stretching the unit cell from a = 1.414 to a = 1.426 Å. Accordingly, the interference between the intraband and interband plasmons strongly increases. A similar interference has been reported in undeformed 5AGNR-arrays doped by positive Fermi energy shifts larger than 0.4 eV [31]. However, such doping values seem to be impractical for current GNR applications.

Conclusion
We have presented a full ab initio modeling, based on groundstate local density calculations, followed by linear response theory, within the RPA, to explore the tunability properties of plasmon excitations in (infinitely periodic) semiconducting (armchair) and semimetallic (zigzag) arrays of GNRs, with ideal symmetric edges, passivated by hydrogen atoms.
All the tested structures are characterized by two interband plasmons at energies larger than 2 eV, which are analogous to the π and π-σ plasmons of graphene. Their peak positions and dispersions are mostly influenced by the GNR width.
At energies smaller than 2 eV, two more intriguing collective excitations appear, which correspond to recently reported edge and surface plasmons [26]. These modes are strongly sensitive not only to the extrinsic conditions, but also to a bunch of geometrical or conformation parameters, such as the width, chirality and unit-cell extension of each GNR, as well as the in-plane vacuum distance between two contiguous GNRs. In ZGNs the absence of a proper LDA band gap, prevents low energy interband excitation from producing a well-defined edge-plasmon structure.
It is worth mentioning that some recent tight-binding models have investigated the role of edge roughness due to asymmetric defects [61,62], which is at present impractical by TDDFT. These studies suggest the edge-plasmon resonances of narrow GNRs are shifted to lower wavelengths and the corresponding plasmon propagation suffers from higher losses with respect to the ideal case, presented here.
We expect that our findings, combined with non-ab initio approaches suitable for the device scale, may open new strategies to construct materials with plasmonic resonances that will be tunable to a specific demand in both the UV-vis and THz regimes, by altering the chemical doping, electronic gating, and also by means of a careful choice of the geometry.